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GLOBAL  DISCRETE  ARTIFICIAL  BOUNDARY  CONDITIONS  FOR 
TIME-DEPENDENT  WAVE  PROPAGATION* * * § 

V.  S.  RYABEN’KIlt,  S.  V.  TSYNKOV§1f,  AND  V.  L  TURCHANINOVH 

Abstract.  We  construct  global  artificial  boundary  conditions  (ABCs)  for  the  numerical  simulation  of 
wave  processes  on  unbounded  domains  using  a  special  non-deteriorating  algorithm  that  has  been  developed 
previously  for  the  long-term  computation  of  wave-radiation  solutions.  The  ABCs  are  obtained  directly  for  the 
discrete  formulation  of  the  problem;  in  so  doing,  neither  a  rational  approximation  of  “non-reflecting  kernels,” 
nor  discretization  of  the  continuous  boundary  conditions  is  required.  The  extent  of  temporal  nonlocality  of 
the  new  ABCs  appears  fixed  and  limited;  in  addition,  the  ABCs  can  handle  artificial  boundaries  of  irregular 
shape  on  regular  grids  with  no  fitting/adaptation  needed  and  no  accuracy  loss  induced. 

The  non-deteriorating  algorithm,  which  is  the  core  of  the  new  ABCs,  is  inherently  three-dimensional,  it 
guarantees  temporally  uniform  grid  convergence  of  the  solution  driven  by  a  continuously  operating  source 
on  arbitrarily  long  time  intervals,  and  provides  unimprovable  linear  computational  complexity  with  respect 
to  the  grid  dimension.  The  algorithm  is  based  on  the  presence  of  lacunae,  i.e.,  aft  fronts  of  the  waves, 
in  wave-type  solutions  in  odd-dimension  spaces.  It  can,  in  fact,  be  built  as  a  modification  on  top  of  any 
consistent  and  stable  finite- difference  scheme,  making  its  grid  convergence  uniform  in  time  and  at  the  same 
time  keeping  the  rate  of  convergence  the  same  as  that  of  the  non-modified  scheme. 

In  the  paper,  we  delineate  the  construction  of  the  global  lacunae-based  ABCs  in  the  framework  of  a 
discretized  w^ave  equation.  The  ABCs  are  obtained  for  the  most  general  formulation  of  the  problem  that 
involves  radiation  of  waves  by  moving  sources  (e.g.,  radiation  of  acoustic  waves  by  a  maneuvering  aircraft). 
We  also  present  systematic  numerical  results  that  corroborate  the  theoretical  design  properties  of  the  ABCs’ 
algorithm. 

Key  words,  artificial  boundary  conditions,  wave  propagation,  lacunae,  non-deteriorating  method 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Numerical  simulation  of  wave  phenomena  on  unbounded  domains  (e.g.,  radiation 
and/or  scattering  of  acoustic  or  electromagnetic  waves),  often  encounters  the  following  well-recognized  dif¬ 
ficulty.  As  no  computer  can  either  handle  infinite  arrays  of  data  or  perform  infinite  numbers  of  arithmetics 
operations,  the  discrete  approximation  of  the  problem  has  to  be  made  finite  (i.e.,  finite-dimensional).  Con¬ 
sequently,  the  original  infinite  domain  has  to  be  truncated  and  special  artificial  boundary  conditions  (ABCs) 
have  to  be  developed  as  a  closure  for  the  resulting  finite  formulation. 

The  literature  on  the  subject  of  ABCs  is  very  extensive,  see,  e.g.,  review  papers  by  Givoli  [1]  and 
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Tsynkov  [2],  as  well  as  another  recent  review  by  Hagstrom  [3],  which  is  geared  more  specifically  toward  wave 
propagation  problems.  In  the  current  study,  we  are  going  to  focus  on  genuinely  unsteady  (as  opposed  to 
time-harmonic)  wave  phenomena  to  be  computed  in  the  time  domain.  For  this  type  of  problems,  most  of 
the  ABCs’  research  to  date  has  been  done  in  the  framework  of  simple  approximate  local  methods  based, 
e.g.,  on  quasi-one-dimensional  characteristics’  decomposition.  These  methods  often  appear  insufficiently 
accurate.  Some  of  the  more  accurate  methods  that  have  recently  gained  attention  are  based  on  the  so-called 
perfectly  matched  layers  (PML),  see  the  original  publications  [4-7]  and  reviews  [8,9].  Unfortunately,  as 
shown  in  [10,11],  these  methods  may  give  rise  to  instabilities  of  a  particular  kind.  The  latter  typically 
manifest  themselves  when  integrating  over  long  time  intervals  and  thus  exacerbate  even  further  the  well- 
known  problem  of  accumulation  of  error  in  long-term  numerical  simulations. 

Among  other  existing  unsteady  ABCs’  approaches,  only  very  few  methodologies  can  guarantee  the  accu¬ 
racy  that  theoretically  would  not  hamper  that  of  the  interior  approximation.  All  of  these  methodologies  are 
non-local,  see,  e.g.,  [12-18],  which  is  characteristic  of  highly-accurate  (ideally,  exact)  ABCs,  The  techniques 
of  this  group  typically  involve  two  “approximating”  steps,  which  are  undertaken  consecutively  when  build¬ 
ing  the  ABCs.  The  first  step  is  the  replacement  of  the  fully  non-local  in  space-time  true  exact  boundary 
conditions,  which  are  most  often  written  using  pseudo-differential  operators,  with  approximate  boundary 
conditions  (still  at  the  continuous  level)  that  would  provide  for  only  a  limited  extent  of  non-locality.  More 
precisely,  this  step  aims  at  limiting  the  temporal  nonlocality  of  the  ABCs,  which  may  be  prohibitively  expen¬ 
sive  in  computations.  This  can  be  achieved,  e.g.,  by  employing  a  rational  approximation  of  the  corresponding 
symbol  (kernel).^  The  first  step  is  then  followed  by  the  second  one,  which  is  the  discretization  of  the  result¬ 
ing  localized  continuous  boundary  conditions.  We  note  that  unless  given  a  special  thorough  attention,  the 
discretization  step  may  lead  to  accuracy  loss  and  again  cause  instability  (this  pertains  to  purely  local  ABCs 
as  well) .  We  also  note  that  all  of  these  techniques  are  restricted  geometrically  to  computational  domains  of 
a  simple  regular  shape,  e.g.,  those  with  spherical  or  planar  boundaries. 

Our  recent  work  [22, 23]  indicates  that  the  issue  of  time-dependent  ABCs  may  be  closely  related  to 
another  problem  that  has  been  mentioned  before  —  the  accumulation  of  numerical  error  during  long  runs. 
This  problem  heis  been  recognized  as  an  outstanding  question  in  computational  PDFs  for  many  years,  since 
the  first  systematic  convergence  studies  for  discrete  approximations  have  been  conducted  in  the  fifties.  Prom 
the  standpoint  of  practical  computing,  deterioration  of  the  solution  over  long  time  intervals  can  be  attributed, 
e.g.,  to  the  mechanism  of  either  numerical  dissipation  or  dispersion  or  both.  Theoretically,  this  phenomenon 
is  conveniently  termed  as  non- uniformity  of  the  grid  convergence  in  time,  and  all  conventional  discrete 
approximations  that  are  currently  in  use  are  known  to  suffer  from  this  deficiency. 

As  our  work  [22,23]  suggests,  the  key  to  resolving  the  question  of  long-term  error  accumulation  may 
be  in  precisely  following  the  physical  nature  of  the  original  problem  when  building  a  numerical  algorithm. 
Namely,  it  is  known  that  waves  in  three  dimensions  have  aft  (or  trailing)  fronts.  This  is  a  manifestation  of 
the  so-called  Huygens’  principle  or  more  generally,  the  presence  of  lacunae  in  wave- type  solutions  in  odd- 
dimension  spaces.  Using  this  property,  we  have  been  able  to  develop  a  modification  to  any  consistent  and 
stable  finite-difference  scheme  that  approximates  a  Cauchy  problem  for  the  wave  equation  so  as  to  make  its 
grid  convergence  uniform  in  time  on  arbitrarily  long  intervals.  The  uniform  temporal  convergence  of  the 
algorithm  has  been  proven  theoretically  along  with  its  optimal  computational  complexity  (i.e,,  linear  with 

^In  fact,  a  wide  variety  of  purely  local  ABCs  (i.e.,  local  in  both  space  and  time)  can  be  obtained  via  rational  approximation 
of  symbols  as  well;  this  approach  has  been  known  for  two  decades,  see  [19-21].  The  general  trend  is  that  the  more  of  the 
nonlocal  nature  of  exact  ABCs  is  compromised,  the  less  accuracy  one  can  expect  from  the  resulting  approximate  methodology. 
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respect  to  the  grid  dimension).  The  rate  of  temporally  uniform  grid  convergence,  see  [22,23],  remains  the 
same  as  that  of  the  original  non-modified  scheme.  These  results  apply  to  the  general  case  of  moving  sources 
of  waves  that  operate  continuously  in  time.  As  an  example,  we  show  in  [23]  that  the  linearized  flow  around 
a  maneuvering  aircraft  can  be  interpreted  in  this  framework. 

At  least  as  important,  the  procedure  of  [22,23]  allows  one  to  replace  the  original  infinite  domain  of  the 
initial- value  problem  by  a  finite  computational  domain  that  would  facilitate  the  construction  of  a  finite¬ 
dimensional  discretization.  As  will  be  seen  from  the  discussion  in  the  current  paper,  the  latter  replacement 
leads  to  obtaining  highly  accurate  non-local  unsteady  ABCs  for  combined  problems  that  may  include  complex 
phenomena  on  a  bounded  interior  domain  but  reduce  to  the  homogeneous  wave  equation  in  the  far  field. 
Similarly  to  the  case  analyzed  in  [22,23]  the  interior  domain  may  be  moving.  The  ABCs  are  built  directly 
for  the  specific  interior  approximation  used  and  in  this  sense  can  be  considered  its  most  natural  extension. 
We  emphasize  that  they  involve  neither  of  the  two  common  approximating  steps  (rational  approximation 
of  the  symbol  followed  by  discretization)  that  have  been  mentioned  before  in  connection  to  some  existing 
ABCs’  methodologies. 

Unlike  in  all  other  methods  available  in  the  literature,  the  extent  of  temporal  nonlocality  of  the  unsteady 
ABCs  based  on  the  technique  [22,23]  is  bounded  naturally  for  all  times  because  of  the  explicit  use  of  lacunae. 
This  bound  is  not  a  consequence  of  any  approximation.  Moreover,  as  these  ABCs  are  obtained  directly  for 
the  specific  finite-difference  scheme,  the  issue  of  discretization  of  the  boundary  conditions,  which  has  been 
shown  to  cause  problems  before,  simply  does  not  arise  in  this  framework.  Besides,  the  new  ABCs  possess 
full  geometric  universality,  i.e.,  can  handle  any  shape  of  the  external  artificial  boundary  with  equal  ease  on 
a  regular  Cartesian  grid  with  no  fitting/adaptation  required  and  no  accuracy  loss  caused. 

The  rest  of  the  paper  is  organized  as  follows.  In  Section  2,  we  provide  for  a  brief  outline  of  the  phe¬ 
nomenon  of  lacunae  in  wave  radiation  solutions,  and  show  how  one  can  use  those  to  obtain  a  non-deteriorating 
algorithm  for  the  long-term  numerical  integration  of  the  corresponding  problems.  In  addition  to  the  theoret¬ 
ical  justification,  we  also  include  in  this  section  several  computational  demonstrations  of  the  properties  of  the 
aforementioned  algorithm.  In  Section  3,  we  describe  in  detail  the  construction  of  the  global  finite-difference 
lacunae-based  ABCs,  and  briefly  comment  on  how  the  proposed  construction  fits  into  the  general  framework 
of  discrete  time-dependent  boundary  conditions  developed  by  Ryaben’kii  in  [24].  Finally,  Section  4  contains 
an  extensive  set  of  numerical  experiments  with  the  new  ABCs  for  the  wave  equation.  The  experiments 
are  conducted  for  finite-difference  schemes  of  different  orders  of  accuracy,  different  laws  of  motion  for  the 
waves’  sources  (uniform,  as  well  as  non-uniform),  and  different  interior  models  that  require  closure  by  the 
homogeneous  wave  equation  in  the  far  field.  These  experiments  corroborate  the  theoretical  design  properties 
of  the  ABCs’  algorithm. 
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2.  Lacunae  and  Non-Deteriorating  Numerical  Integration. 

2.1.  Lacunae  of  the  Wave  Equation.  We  consider  a  Cauchy  (initial- value)  problem  for  the  three- 


dimensional  wave  equation,  x  =  {xi^x^-^xz)- 


df- 


-  c 


dx\  dx\  dx\ )  ’  ’ 


(2.1) 


t-0 


d(f 

dt 


=  0  . 
t=o 


(2.2) 


(The  limitation  of  having  homogeneous  initial  conditions  (2.2)  can  be  alleviated,  see  [22,23].)  The  problem 
(2.1),  (2.2)  is  solved  on  the  domain  S{t)  C  which  has  finite  diameter  d  for  all  times  t  >  0;  other  than  that 
the  domain  S{t)  may  travel  in  space  according  to  an  arbitrary  law  of  motion  except  that  its  maximum  speed 
k  is  required  to  be  “subsonic:”  k  <  c.  The  solution  (p(x,t)  is  driven  by  the  continuously  operating  source 
f{x,t),  f{x,0)  =  0,  and  we  require  that  \/t  >  0:  supp/(a!;,t)  C  S{t).  In  other  words,  we  study  the  radiation 
of  waves  by  a  source,  which  is  compactly  supported  in  space  for  all  times.  The  solution  is  of  interest  for  us 
also  on  a  compact  domain,  which  we  call  S{t)]  it  fully  contains  the  source  and  follows  its  motion  if  there  is 
motion.  This  is  a  simplified  model  for  many  interesting  physical  phenomena  that  are  more  complex  in  their 
nature.  As  we  shall  see,  this  model  is  very  useful  when  constructing  the  ABCs  for  a  variety  of  problems. 

For  every  (x,  t),  the  solution  (p  =  (p(x,t)  of  problem  (2.1),  (2.2)  is  given  by  the  Kirchhoflf  integral: 


(fix^t) 


g<ct 


Q 


(2.3) 


where  |  =  (^i,^2,^3),  ^  =  \/{^i  ”  6)^  +  (^2  -  6)^  +  =  d^id^2d^z.  If  we 

assume  for  a  moment  that  the  right-hand  side  (RHS)  f{x^t)  is  compactly  supported  in  space-time  on  the 
domain  Q  x  [0,-hoo),  then  formula  (2.3)  immediately  implies  that 


(p{x,t)=0  for  {x,t)  e  Pi  {(ic,t)||a;  -  <  c(i  -  ^),  t  >  ^}  .  (2.4) 

imeQ 

The  region  of  space-time  defined  by  formula  (2.4)  is  called  lacuna  of  the  solution  (p  =  (p{x,t).  This  region 
is  obviously  obtained  as  the  intersection  of  characteristic  cones  of  equation  (2.1)  once  the  vertex  of  the  cone 
sweeps  the  support  of  the  RHS:  supp  /  C  Q.  Prom  the  standpoint  of  physics,  the  lacuna  represents  that  part 
of  space-time  where  the  waves  generated  by  the  sources  /,  supp  /  C  (J,  have  already  passed  and  the  solution 
has  become  zero  again.  (Sometimes,  the  name  “secondary  lacuna”  is  used  to  distinguish  it  from  the  primary 
lacuna,  which  is  the  area  where  the  waves  have  not  reached  yet.)  The  phenomenon  of  lacunae  is  inherently 
three-dimensional.  The  interior  surface  of  the  lacuna  represents  the  trajectory  of  aft  (trailing)  fronts  of  the 
waves.  The  presence  of  aft  fronts  in  odd-dimension  spaces  is  known  as  the  Huygens’  principle,  as  opposed 
to  the  so-called  wave  diffusion  which  takes  place  in  even-dimension  spaces. 

2.2.  Computation  with  a  Compactly  Supported  Source.  Assume  now  that  the  moment  of  in¬ 
ception  of  the  source  f{x,t)  is  to  (in  particular  it  may  be  to  =  0);  at  this  moment  the  domain  5(to)  of  the 
RHS  /(x,  t)  occupies  a  position  in  space  which  is  schematically  represented  by  the  interval  [^1,^2]  of  size  d 
on  Figure  2.1.^  Assume  also  that  by  the  time  h  >  to  this  source  ceases  to  operate,  which  makes  the  RHS  of 

^Throughout  this  section  we  will  be  using  schematic  one-dimensional  illustrations  always  keeping  in  mind,  however,  that 
the  eictual  problem  is  three-dimensional. 
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Fig.  2.1.  ID  schematic  representation  for  the  fronts  of  waves  generated  by  a  compactly  supported  source. 

equation  (2.1)  compactly  supported  in  both  space  and  time:  supp/  C  Q  =  a:  €  S{t),  Iq  <t  <  ti}. 

Clearly,  by  the  time  h  the  domain  S(ti)  can  be  displaced  from  its  initial  location  no  further  than  the  distance 
^(^1  “^o)  is  each  direction,  which  is  schematically  represented  on  Figure  2.1  by  the  boundaries  of  the  interval 
[Bi,B2]  of  size  d  +  2k{ti  -  ^o)-  Starting  from  ^  no  new  waves  will  be  generated,  and  those  generated 
prior  to  ti  will  continue  traveling  in  space  and  thus  will  eventually  leave  the  domain  S (t)  completely,  because 
their  speed  of  propagation  c  is  higher  than  that  of  the  domain,  k.  The  moment  t2  when  this  happens,  i.e., 
when  the  solution  again  becomes  zero  on  5(t),  is  easy  to  calculate,  see  Figure  2.1.  By  this  moment, 

the  domain  5(^2)  travel  no  further  than  the  interval  [Ci,C2]  of  size  d  +  2k{t2  ~  to)^  and  we  need  to 
assume  that  the  aft  fronts  will  also  be  exactly  at  the  boundaries  of  this  interval  at  t  =  ^2,  which  immediately 
yields 

ip{x,t)  =  0,  for  xeS{t),  t>t2=to  + - ^  ^  _  (2.5) 

Estimate  (2.5)  is  fundamental.  It  essentially  says  that  once  we  need  to  calculate  the  solution  (p(x,  t)  on  S{t) 
and  the  sources  are  compactly  supported  in  space-time:  supp  /  C  Q  {{x,t)\x  e  S{t),  to  <t  <ti}y  then  we 
may  stop  the  calculation  aX  t  =  t2  because  afterwards  the  solution  on  S{t)  will  be  zero  anyway.  This  means, 
in  particular,  that  if  the  solution  is  calculated  using  a  discrete  method,  e.g.,  a  finite-difference  scheme,  then 
no  new  error  will  be  accumulated  after  t  =  ^2*  The  constants  in  both  consistency  and  stability  estimates 
of  the  scheme  (see  [22,23]  and  below  for  detail)  will  now  depend  on  the  time  interval  Tint  = 
rather  than  final  time  Tfinai,  which  for  the  case  of  a  compactly  supported  RHS  simply  becomes  immaterial. 

Besides,  once  we  stop  the  calculation  at  ^2  ==  +  Tint?  we  realize  that  during  the  time  interval  Tint  that 

has  passed  since  the  beginning  t  ~  to^  xio  waves  could  have  traveled  in  space  further  than  the  boundaries 
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of  the  interval  [I>i,D2]  of  size  d  +  2cTint,  see  Figure  2.1.  Beyond  this  region  the  solution  is  zero  because 
this  is  the  area  of  the  primary  lacuna.  Therefore,  even  though  the  original  problem  was  formulated  on  an 
infinite  domain,  we  can,  in  fact,  calculate  the  solution  on  a  finite  domain  of  size  d  +  2cTint  with  zero 

external  boundary  conditions  (of  the  Dirichlet  type). 

The  transition  from  the  infinite  domain  to  a  finite  one  does  not,  obviously,  come  “at  no  charge.”  One 
can  rather  say  that  it  comes  at  the  expense  of  having  the  computational  domain  [Di ,  D2]  larger  than  the 
actual  domain  of  interest  S{t).  However,  the  size  of  the  “redundant”  portion  of  [£^1,1)2]  can  be  further 
reduced  by  observing  that  all  we  have  to  do  is  make  sure  that  by  i  by  the  moment  the  last  waves 

generated  by  f{x,t),  supp/  C  Q,  leave  5(t),  no  new  waves  can  enter  S(t).  This  can  be  guaranteed  either 
as  indicated  previously,  by  placing  the  outer  boundary  sufficiently  far  away  so  that  no  waves  from  /(x,t), 
supp  /  C  Q,  can  even  reach  it  until  t  =  t2  =  to  -\-  Tint  or  alternatively,  by  placing  it  closer  and  thus  allowing 
for  reflections,  but  still  not  too  close  so  that  no  reflected  waves  can  come  back,  i.e.,  reach  5(t),  by  ^  =  ^2- 
The  size  of  the  new,  smaller,  computational  domain  [Ei^E2]  with  reflecting  outer  boundary,  see  Figure  2.1, 
can  be  estimated  easily.  The  minimum  size  Z,  see  Figure  2.1,  is  found  by  requiring  that  the  reflected  waves, 
which  travel  with  the  same  speed  c  but  in  the  opposite  direction,  reach  the  boundary  of  [C'i,C'2],  i.e.,  the 
utmost  possible  location  of  5(^2),  by  the  exact  same  moment  of  time  t  =  t2  when  the  aft  fronts  leave  S{t). 
This  immediately  yields 

Z  =  d+(^•  +  c)Tint  . 

By  comparing  the  value  of  Z  from  (2.6)  with  the  size  of  [Di,T)2]j  which  is  d  -f-  2cTint5  we  conclude  that  the 
extra  size  of  the  computational  domain  beyond  d  can  be  reduced  by  up  to  a  factor  of  2  (when  k  —  0)  in  each 
coordinate  direction. 

We  also  note  that  in  fact  any  well-posed  boundary  condition  can  be  specified  at  the  reflecting  outer 
boundary  of  [FJi,  E2].  The  reason  is  that  this  boundary  is  intentionally  positioned  so  that  the  reflections  are 
not  going  to  have  any  effect  on  the  solution  inside  S{t)  anyway.  A  particularly  convenient  way  to  treat  the 
boundary  of  ,  E2]  will  be  to  set  the  periodic  boundary  conditions  there.  In  so  doing  the  three-dimensional 
rectangular  domain  becomes  a  three-dimensional  toroidal  surface  (the  opposite  faces  of  the  rectangle  are 
identified  with  one  another)  and  we  only  have  to  keep  in  mind  that  the  reflected  waves  will  now  need  to 
be  interpreted  as  those  that  leave  the  domain  on  one  side  and  enter  it  from  the  opposite  side.  This  new 
interpretation  obviously  brings  no  change  to  the  foregoing  considerations  that  led  to  the  size  estimate  (2.6). 
However,  for  the  case  of  a  continuously  operating  traveling  source  that  we  analyze  below,  periodicity  implies 
that  the  motion  of  the  source  can  also  be  formally  considered  on  the  toroidal  surface,  which  makes  the 
computational  setup  much  simpler. 

2.3.  Computation  with  a  Continuously  Operating  Source.  Both  foregoing  observations  finite 
time  interval  Tint  and  finite  spatial  domain  [Ei ,  £^2]  needed  for  calculating  the  solution  driven  by  the  sources 
f{x,  t):  supp  f  CQ  =  {{x,t)\x  e  S{t),  to  <t  <ti},on  the  domain  S{t)  —  are  crucial  for  the  original  case 
of  a  continuously  operating  source  /(x,t):  supp/  C  {{x,t)\x  e  S{t)^  t>0}.  In  this  case  we  first  take  a 
parameter  T  >  0  and  introduce  a  smooth  even  compactly  supported  function  ©(t),  t  E  M,  of  a  “hat”  type: 

©(t)  =  0  ,  lt|  >  T  , 
e{t)  =  e{~t) , 

©(t)  =  1  ,  t  e  [-(jT,  crT]  ,  0  <  a  <  1  , 
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which  obviously  generates  a  partition  of  unity: 

oo 

l  =  ^0{t-(l  +  a)Ti)  ,  ^>0, 

with  the  overlap  size  (1  —  (t)T.  Then,  we  represent  the  right-hand  side  f{x,t)  of  equation  (2.1)  in  the  form 

CXD  oo  CO 

fix, t)  =fix,t)^&it-il  +  a)Tj)  =  ^ ©((  -  (1  +  a)Tj)fix,t)  =  fj{x,t, T)  .  (2.7) 

j=0  j=0  j=o 

Clearly,  for  each  fj{x,t,T)  =  G{t  ~  (1  +  a)Tj)f{x,t)y  j  —  0, 1, . . . ,  we  have 

snpp  f j{x,t,T)  C  {{x,t)\x  e  S{t),  {j -l  +  crj)T  <t  <{j -\-l-\-crj)T}  .  (2.8) 

Due  to  the  linear  superposition,  the  overall  solution  (p{x^t)  of  equation  (2.1)  will  be  given  by  the  sum  of 
individual  contributions  from  /j(a:,  T),  j  =  0, 1, . . . : 

oo 

^ix,t)  =  '£<Pj{x,t,T)  ,  (2.9) 

J=0 


where  each  contribution  (pj{Xyt,T)  solves  the  following  sub-problem: 


_  2  { 

dt^  ^  \  dxl  dx2  dx\  ) 


fj{Xj t, T)  , 


h=0'-i+<7j)T 


=  0,  j  =  0,l,2,... 


(2.10) 


Notice  that  each  (pj(x,t,T),  j  =  0, 1, . . . ,  can  be  calculated  absolutely  independently  of  the  others  and  that 
the  corresponding  source  term  //(x,  t,  T)  is  a  function  compactly  supported  in  both  space  and  time,  see  (2.8). 
Consequently,  according  to  (2.5),  if  we  interpret  to  and  h  as  ( j  -  1  +  (7j)T  and  ( j  +  1  +  (rj)T,  respectively 
(see  Figure  2.1),  then  we  can  conclude  that  every  ipj{x,t^T)  of  (2.9)  needs  to  be  calculated  only  during  a 
finite  interval  of  time  Tint  =  ^  is  important  to  realize  that  this  interval  does  not  depend  on  the 

actual  moment  of  time  t. 

Moreover,  even  so  the  series  (2.9)  is  formally  infinite,  it  is  easy  to  see  that  for  any  t  >  0,  x  E  S{t)j  it 
contains  only  a  finite  fixed  number  of  non-zero  terms.  First  of  all,  because  of  the  causality,  ipj{x^t,T)  =  0 
for  x  E  S{t)  if  t  <  (j  -  1  +  aj)T.  In  other  words,  for  a  given  moment  of  time  t,  the  contribution  of  all  those 
/j(x,t,T)  that  are  active  only  at  subsequent  moments  of  time,  is  obviously  zero.  A  somewhat  less  trivial 
observation  is  that  because  of  the  lacunae  the  contribution  of  the  ‘^sufficiently  retarded”  terms  fj{x,t,T)  to 
the  overall  solution  at  a  given  time  level  t  will  be  zero  as  well.  More  precisely,  ipj{x^t,T)  =  0  for  x  E  S{t) 
if  (i  -  1  +  (rj)T  <  t  -  Tint.  This  follows  immediately  from  (2.5)  assuming  that  to  -  {j  -  I  +  (^j)T  and 
=  (j  +  1  +  aj)T.  Consequently,  instead  of  (2.9)  we  can  write: 


P2 

ip{x,t)  =  '£^jix,t,T),  xeS{t),  (2.11) 

j=Pl 

where  pi  =  +  l)] ,  P2  =  (f  +  l)]  ’  integer  part.  The  expressions  for 

Pi  and  p2  indicate  that  we  will  always  have  either  pi  —P2-  [(i+^yr]  Pi  =  P2  -  [  (i+aV]  “  Therefore, 
the  number  of  terms  p  =  p2  —  {pi  —  V)  in  the  sum  (2.11)  will  never  exceed  +  2.  As  Tint  does  not 
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depend  on  t  we  conclude  that  neither  does  the  foregoing  upper  bound  for  p.  As  such,  the  number  of  terms 
in  the  sum  (2.11)  can  always  be  considered  finite  and  fixed.  Altogether  we  obtain  that  Vt  >  0  the  solution 
ip{x,t),  X  e  5(t),  is  composed  of  a  finite  non-increasing  number  p  or  additive  terms,  and  each  of  the  latter 
needs  to  be  taken  into  account  only  during  a  finite  non-increasing  interval  of  time  Tint- 

In  the  perspective  of  numerical  computation,  the  latter  consideration  translates  into  temporally  uniform 
grid  convergence  of  the  discrete  algorithm.  Indeed,  assume  that  we  are  integrating  equation  (2.1)  by  means 
of  a  finite-difference  scheme  with  the  order  of  accuracy  0{h^)j  where  h  is  a  general  notion  for  the  grid  size 
and  a  >  0.  Then,  the  discrete  solution  (p^^\x^t)  converges  to  the  continuous  solution  as  the  grid 

size  decreases: 

||^W(a:,t)  - ^(x,t)||  <K-h^  ,  t&  [O.Tfinai]  ,  (2-12) 

here  Tfinai  is  the  total  integration  time.  Inequality  (2.12)  is  a  generic  convergence  estimate;  it  holds  provided 
that  the  RHS  /(x,  t)  of  equation  (2.1)  is  sufficiently  smooth.  A  detailed  discussion  on  the  smoothness  require¬ 
ments  for  f{x,t)  can  be  found  in  [23],  along  with  the  specific  consistency /stability /convergence  estimates  in 
the  norms  that  would  take  into  account  a  particular  smoothness  level. ^ 

The  constant  K  in  inequality  (2.12)  does  not  depend  on  the  grid.  It  is,  however,  known  to  depend  on 
the  actual  RHS  f{x,t),  as  well  as  the  final  time  Tfinai*  AT  =  /^(Z,  Tfinai)-  The  dependency  of  K  on  Tjfinai  is 
typically  a  grow^th,  and  sometimes  this  growth  may  be  rapid.  This  means  that  even  so  on  any  fixed  interval 
[0,  Tfinai]  the  scheme  converges  as  h  — >  0,  to  obtain  the  same  level  of  accuracy  on  a  larger  [0,  Tfinai]  one  may 
need  to  take  a  finer  overall  grid  ahead  of  time.  Thus,  the  convergence  appears  temporally  non-uniform.  On 
the  language  of  practical  computing,  this  phenomenon  can  be  interpreted  as  the  accumulation  of  numerical 
error  over  long  runs.  This  issue  has  been  long  acknowledged  unresolved  in  the  literature. 

The  situation  changes  dramatically  if,  instead  of  the  straightforward  time-dependent  integration,  we 
first  use  the  foregoing  lacunae-based  representation  (2.11)  of  the  solution  ip{Xyt).  In  so  doing,  for  each 
j  =  pi,...  ,_p2,  we  still  integrate  the  corresponding  sub-problem  (2.10)  using  the  same  finite-difference 
scheme  as  before.  However,  convergence  estimate  for  the  scheme  then  becomes: 

)  {X,  t,  T)  -  ^j{x,  t,  T)  I  <  Kj  •  /»“  , 

X  e  S{t)  ,  t  e[{j  -l-\-  (Tj)T,  (j  -  1  -h  aj)T  -h  Tint]  • 

A  very  important  circumstance  is  that  unlike  K  in  estimate  (2.12),  the  constant  Kj  in  (2.13)  for  each  j 
depends  on  Tint  rather  than  Tfinai:  Kj  =  Kj{fj,Tint).  Keeping  in  mind  that  each  ipf\x,t,T)  can  be 
computed  independently  of  the  others,  and  using  linear  superposition  (formula  (2.10)),  we  then  easily  obtain 
instead  of  (2.12): 

||(^W(x,t) -(p(x,i)||  ,  x£S{t),t>0,  (2.14) 

where  K  =  A"(/,Tint).  Note  that  p  is  fixed  and  does  not  increase  with  t,  and  K  now  depends  on  Tint  rather 
than  Tfinai,  where  Tint  is  also  fixed  and  does  not  increase  with  t  Therefore,  estimate  (2.14)  implies  temporally 
uniform  grid  convergence  of  the  discrete  lacunae-based  algorithm  on  arbitrarily  long  time  intervals  or  in  other 
words,  for  t  >  0.  A  detailed  formal  proof  of  this  result  that,  again,  involves  specific  norms,  can  be  found 
in  [23].  Here  we  only  need  to  add  that  for  each  inequality  (2.13)  to  hold,  the  corresponding  /j(aj,t,T), 

^Smoothness  of  the  source  terms  will  also  be  important  when  constructing  the  lacunae-based  ABCs,  see  Sections  3  and  4. 


8 


j  =  pi,. ..  ,p2>  lias  to  possess  the  same  regularity  as  that  required  for  the  original  scheme  to  converge.  This 
explains  why  choosing  the  partition  (2.7)  smooth  with  overlaps  was  very  important. 

From  the  standpoint  of  practical  computing,  temporally  uniform  grid  convergence  implies  that  the 
numerical  error  will  not  get  accumulated  beyond  some  predetermined  bound  for  as  long  as  the  computation 
needs  to  be  performed;  and  once  the  grid  is  refined  the  aforementioned  bound  will  also  drop  in  accordance  with 
the  specific  rate  0{h^).  This  is  clear  because  when  the  calculation  is  stopped  for  a  given  term  (p^^\x,t,T) 
after  the  interval  Tint  has  elapsed,  the  error  will  not  be  accumulated  any  further,  and  the  number  of  terms  p 
that  need  to  be  taken  into  account  is  fixed  and  non-increasing.  Thus,  we  have  obtained  a  non-deteriorating 
numerical  algorithm  for  integration  of  the  wave  equation  over  arbitrarily  long  times.  Let  us  emphasize  that 
it  can  be  built  as  a  modification  of  any  consistent  and  stable  finite-difference  scheme,  and  that  it  preserves 
the  original  rate  of  convergence  of  the  scheme  while  making  the  convergence  uniform  in  time. 

Besides,  let  us  assume,  for  example,  that  the  original  finite-difference  scheme  has  linear  computational 
complexity  with  respect  to  the  grid  dimension,  which  is  typical  for  explicit  schemes.  Then,  it  is  easy  to  see 
that  the  modified  lacunae-based  algorithm  will  also  have  linear  computational  complexity  with  respect  to  the 
grid  dimension.  Indeed,  this  immediately  follows  from  the  fact  that  each  term  (x,  t,  T)  is  computed  using 
the  original  scheme  on  a  compact  domain  of  size  Z  (see  Figure  2.1)  during  a  finite  fixed  interval  of  time  Tint, 
and  the  number  of  terms  p  is,  again,  fixed  and  non-increasing.  We  should  note  that  for  the  type  of  problems 
that  we  are  studying  linear  complexity  with  respect  to  the  grid  is,  in  fact,  optimal,  i.e.,  unimprovable. 

2.4.  Computation  Using  Continuous  Time  Marching.  The  following,  and  last,  step  in  building 
the  lacunae-based  algorithm  for  long-term  numerical  integration  of  the  wave  equation  is  to  realize  that  for 
implementing  formula  (2.11)  we  do  not  necessarily  need  to  compute  each  term  (pj{x,t,T)  independently  of 
the  others.  Instead,  we  can  implement  the  algorithm  in  a  way  similar  to  the  standard  time-marching  by 
means  of  a  finite-difference  scheme.  For  that,  we  will  need  to  use  the  aforementioned  periodic  boundary 
conditions  on  the  outer  boundaries  of  the  auxiliary  domain  [Ei,E2]j  see  end  of  Section  2,2  and  Figure  2.1. 

The  first  key  observation  that  we  make  here  is  that  once  the  motion  of  the  wave  sources,  as  well  as 
the  propagation  of  waves  themselves,  are  considered  on  a  three-dimensional  toroidal  surface,  rather  than  on 
the  genuine  M^,  then  for  every  portion  of  the  RHS  fj(x,t,T)  it  does  not  really  matter  where  on  the  period 
this  source  is  located,  or  where  it  starts  its  motion  from,  at  to  =  ( j  “  1  +  does  not  have  to  be 

exactly  “in  the  middle”  as  shown  on  Figure  2.1,  because  all  locations  on  the  period  (i.e.,  toroidal  surface) 
are  equivalent.  All  we  have  to  worry  about  is  that  by  the  time  t2  =  (i  “  1  +  <^i)^  +  ^int  the  waves  generated 
by  fj{x,t,T),  see  (2.8),  will  have  left  the  domain  5(f),  and  that  no  waves  could  have  re-entered  this  domain 
during  [fo,f2]-  And  this  will  be  exactly  the  case  because  the  size  Z  =  d  +  {k  c)Tint,  see  (2.6),  has  been 
chosen  sufficiently  large  to  provide  for  that.  Since  we  always  assume  (for  simplicity)  that  the  period  Z  is 
the  same  in  all  coordinate  directions,  then  we  only  need  to  formally  consider  /j(x,f,T)  instead  of  fj{x,t,T) 
and  accordingly,  ipj{x,t,T)  instead  of  (pj{x,t,T),  where  x  =  {xi,X2,xz),  and  Xi  =  Xi  —  [^]  Z,  i  =  1,2,3. 

Next,  we  shall  analyze  formula  (2.11)  from  a  slightly  different  point  of  view.  To  begin  with,  we  notice 
that  on  the  initial  stage  of  computation,  i.e.,  when  t  is  small,  the  lower  summation  limit  pi  may  turn  out 
negative.  Basically,  it  does  not  create  any  inconsistency  and  does  not  cause  any  problem  because  f{x,t)  —  0 
for  f  <  0  anyway.  In  fact,  we  can  simply  disregard  all  negative  j’s  in  the  sum  (2.11)  for  small  f’s  and 
initially  consider  the  summation  instead  of  (2.11).  “Initially”  here  means  till  the  actual 

expression  pi  =  +  l)]  becomes  positive.  It  is  easy  to  see  that  the  computation  on  this  initial 

stage  is  equivalent  to  the  conventional  time-marching  of  the  wave  equation  (2.1)  on  the  domain  [Ei^Ez]  of 
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size  Z  (see  Figure  2.1)  with  periodic  boundary  conditions.  Indeed,  all  we  do  here  is  simply  take  into  account 
one  component  of  the  source  /j  after  another.  Due  to  linear  superposition,  this  amounts  to  the  continuous 
integration  of  the  wave  equation  driven  by  /  =  /j  from  t  —  0  till  the  actual  time  t.  We  also  note  that 
the  duration  of  the  initial  stage  is,  obviously,  Tint-  And  the  period  Z,  see  (2.6),  has  been  chosen  sufficiently 
large  so  that  for  the  time  interval  of  length  Tint  there  will  be  no  difference  on  the  domain  S  (t)  between  the 
solution  (p{x,t)  computed  in  the  periodic  setting  and  the  solution  that  one  could  have  possibly  computed 
with  no  periodization  (see  Section  2.2). 

As  soon  as  the  time  interval  Tint  =  has  elapsed  since  the  inception  moment  t  =  0,  the 

computation  enters  its  regular  (as  opposed  to  initial)  stage.  This  regular  stage,  which  can,  in  fact,  be 
continued  for  as  long  as  necessary,  is  characterized  by  the  positive  values  of  pi  (the  first  positive  value  is 
obviously  pi  =  1)  and  finite  non-growing  number  p  =  P2  -  {pi  -  1)  of  terms  in  the  sum  (2.11). 

On  the  regular  stage  of  the  algorithm,  we  continue  marching  equation  (2.1)  with  periodic  boundary 
conditions  in  space.  Obviously,  as  the  time  t  elapses  both  pi  and  p2  in  formula  (2.11)  increase.  The  increase 
of  Pi  and  p2  is  almost  synchronous.  Namely,  as  soon  as  t  reaches  the  value  ( j  —  1  +  (Tj)T  for  a  particular 
integer  j,  a  new  term  (pj  gets  included  into  the  sum  (2.11),  i.e.,  the  upper  summation  bound  p2  changes 
from  its  previous  value  j  —  1  to  the  new  value  j.  Similarly,  as  soon  as  t  reaches  the  value  {j  —  l-\-CFj)T  Tint 
for  a  given  j,  the  term  (pj  drops  from  the  sum  (2.11),  i.e.,  the  lower  summation  bound  changes  from  its 
previous  value  j  to  the  new  value  j  +  1.  As  has  been  mentioned  in  Section  2.3,  in  so  doing  the  variation 
of  the  difference  between  p2  and  pi  never  exceeds  one.  Moreover,  the  temporal  interval  that  precedes  the 
actual  moment  t  and  that  is  taken  into  consideration  by  formula  (2.11)  is  again  Tint-  Consequently,  we  can 
still  compute  everything  in  the  periodic  framework,  because  the  period  Z  (see  (2.6))  is  sufficiently  large  to 
accommodate  the  extent  of  retardation  Tjnt,  and  as  has  also  been  mentioned  it  does  not  matter  where  on 
the  period  the  computation  of  every  given  term  starts. 

Prom  the  standpoint  of  implementation,  when  the  upper  bound  p2  increases  by  one  at  t  =  (j  -  1  4-  crj)T 
nothing  special  needs  to  be  done.  If  we  simply  continue  marching  equation  (2.1)  in  the  aforementioned 
periodic  framework,  then  we  will  automatically  start  taking  into  account  the  new  component  of  the  RHS  fj 
after  t  =  ( j  —  1  +  <tj)T.  The  situation  with  the  lower  bound  pi  is  somewhat  different.  Once  it  has  increased 
by  one  (from  j  to  j  ~\-l)  at  t  =  ( j  —  1  +  o‘j)T  -I-  Tint,  the  term  (pj  no  longer  needs  to  be  included  into  the  sum 
(2.11).  However,  in  contradistinction  to  the  case  of  Section  2.3  when  all  (pj  were  supposed  to  be  computed 
independently  of  one  another,  here  we  cannot  just  stop  the  computation  of  a  given  (pj  at  t  =  {j—l-\-aj)T +Tint 
and  subsequently  say  that  (pj («,  t, T)  =  0  for  x  G  S{t)  and  for  t>  {j- 1  +  <jj)T +Tint-  Indeed,  time-marching 
of  equation  (2.1)  implies  that  all  fragments  of  the  solution  (pj{x,t,T)  are  calculated  together  as  a  sum  and 
cannot  be  explicitly  told  apart.  On  the  other  hand,  if  we  do  nothing  at  t  =  {j  —  l-\-  crj)T  +  Tint  and  continue 
with  the  time-marching,  i.e.,  if  we  do  not  discontinue  the  computation  of  (pj{x,  t,  T)  at  t  =  {j  ~l-\-aj)T +Tint 
and  leave  this  term  in  the  solution  (^(a?,  i),  then  right  after  this  moment  of  time  the  first  waves  generated  by 
fj  at  t  =  {j  —  1  +  crj)T  will  start  re-entering  the  domain  S{t)  having  traveled  all  the  way  across  the  auxiliary 
domain  [Ei,E2]-  In  other  words,  in  the  framework  of  the  continuous  time-marching  with  periodic  boundary 
conditions,  the  term  ipj{x,t,T)  cannot  be  left  in  the  solution  as  it  will  “contaminate”  the  results  on  S{t). 

To  avoid  the  aforementioned  contamination,  i.e.,  to  prevent  the  re-entry  of  waves  into  5(t),  each  term 
(pj{x,  t, T)  needs  to  be  eliminated  from  the  overall  solution  on  the  auxiliary  domain  [FJi,  E2]  when  the  extent 
of  its  retardation  (counted  from  inception)  becomes  exactly  Tint-  For  a  given  j  the  proper  moment  of  time  for 
elimination  of  (pj (x,  t,  T)  is  t  =  (j  —  1  +  crj)T  +  Tint-  Once  we  take  out  (pj {x,  t, T)  at  t  =  (j  —  1  +  crj)T  +  Tint, 
this  term  may  obviously  be  considered  zero  everywhere  on  [Ei,E2]  for  all  subsequent  moments  as  well.  To 
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take  out  the  term  (pj{x,t,T)  we  need  to  interrupt  the  time-marching  at  i  =  (i  -  1  +  crj)T  -h  Tint,  then  go 
back  to  the  inception  moment  of  fj{x,  t, T),  which  is  t  =  (i  -  IH- aj)T,  and  independently  integrate  problem 
(2.10)  for  a  particular  j  on  [£^1,^2]  from  to  —  (i  —  1  +  ^2  =  (i  “  1  +  +  ^int-  The  result  should 

then  be  subtracted  from  the  time-marching  solution  at  t  =  t2  in  the  correct  sense,  i.e.,  both  (p  and  (rather 
their  discrete  counterparts)  should  be  affected.  Alternatively,  we  may  notice  that  when  integrating  problem 
(2.10)  from  to  =  (j  -  1  +  (Tj)T  till  t2  =  ( j  -  1  +  (Tj)T  -h  Tint,  the  wave  equation  will,  in  fact,  be  homogeneous 
on  a  substantial  portion  of  this  time  interval  because  fj{x,  t^T)  :=  0  for  t  >  ti  =  (j  -1-  ,  Consequently, 

instead  of  marching  equation  (2,10)  over  the  entire  time  interval  of  length  Tint,  we  may  actually  march  it 
only  from  to  =  (j  - 1  +  (Tj)T  till  ti  =  {j  -h  l-\-aj)T^  then  Fourier  transform  the  discrete  solution  and  advance 
it  till  ^2  —  (j  -  1  +  ^i)T'  +  Tint  by  raising  the  corresponding  amplification  factors  to  the  appropriate  power. 
Numerically,  this  approach  appears  much  cheaper,  especially  if  it  relies  on  highly  efficient  FFT  subroutines. 

The  new  version  of  the  lacunae-based  algorithm  has  obviously  been  designed  so  that  to  exactly  repro¬ 
duce  the  solution  obtained  with  the  original  version  of  Section  2.3.  The  only  difference  is  in  the  method 
of  computation:  Continuous  time-marching  in  the  periodic  setup  with  cyclic  subtractions  of  the  retarded 
contributions  versus  separate  computation  of  partial  solutions  driven  by  different  components  of  the  RHS. 
Consequently,  the  new  version  will  posses  the  same  properties  as  the  original  one.  Foremost,  it  will  provide 
for  the  temporally  uniform  grid  convergence.  Besides,  it  will  obviously  have  linear  computational  complex¬ 
ity  with  respect  to  the  grid  dimension.  (The  cost  of  the  FFT-based  evolution  in  time  distributed  over  the 
corresponding  number  of  time  steps  is  even  less  than  linear  if  calculated  per  time  step.)  Finally,  the  algo¬ 
rithm  will  be  universal  in  the  sense  that  one  will  be  able  to  build  it  as  a  modification  of  any  consistent  and 
stable  finite-difference  scheme.  It  will  preserve  the  convergence  rate  of  the  original  scheme  while  making  the 
convergence  uniform  in  time. 


2.5.  Numerical  Demonstrations.  To  actually  demonstrate  that  the  lacunae-based  algorithm  is  an 
appropriate  procedure  that  does  deliver  according  to  its  theoretical  design  properties,  we  present  some 
numerical  results  for  the  wave  equation.  For  our  simulations,  we  assume  axial  symmetry  and  employ  the 
(r,  z)  cylindrical  coordinates  so  that  to  account  for  the  three-dimensional  effects  using  two-dimensional 


geometry.  Accordingly,  equation  (2.1)  becomes: 


W 


f{r,z,t)  ,  t>0. 


(2.15) 


The  solution  p  of  equation  (2.15),  as  well  as  the  RHS  /,  are  functions  of  r,  z,  and  t.  The  initial  conditions 


for  equation  (2.15)  remain  homogeneous  as  before,  see  (2.2). 

We  introduce  the  rectangular  auxiliary  domain  [0,R]  x  [— Z/2,  Z/2]  of  variables  (r,  z),  this  domain  is  a 
specific  realization  of  [Ei,E2]  shown  in  Figure  2.1.  The  boundary  conditions  are  periodic  with  the  period  Z 


in  the  z  direction,  and  zero  Dirichlet  sX  r  —  R: 


(p{r,z±  Z,t)  =  (p{r,z,t)  , 

(p{RyZ,t)  =  0  , 


(2.16) 


The  mathematical  formulation  of  the  problem  obviously  requires  no  boundary  conditions  at  r  =  0.  However, 
for  the  purpose  of  subsequently  building  a  discrete  scheme  (see  below)  we  notice  that  the  natural  assumption 
of  cp{r,z^t)  being  a  bounded  smooth  function,  along  with  the  axial  symmetry,  immediately  imply  that 
1^1  =0.  Consequently,  the  Taylor  expansion  for  (p  near  r  =  0  yields: 

cfr  |^_Q 


1  d‘^(p 

p>{r,  ■)  =  v?(0,  ■)  +  2  ^ 


■  -h  0{r^)  , 


r=0 
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which  means  that 


dip  _  d^if 
dr  dr‘^ 


r=0 


.r  +  0(r2)  . 


Substituting  the  latter  expression  into  (2.15)  and  considering  the  limit  r  0,  we  obtain  that  on  the  2;-axiSj 
i.e.,  at  r  =  0,  equation  (2.15)  reduces  to: 


di^if  9 4.\  n  + 


0  . 


(2.17) 


To  assess  the  quality  of  our  numerical  method  we  need  to  build  a  reference  exact  solution  of  problem 
(2.15),  (2.2).  This  solution  is  obtained  using  the  Lorentz’  transform: 

1  ,  kjc  z 

(2.18) 


C  =  - 


yrrPT? 

k/c 


■t- 


sjl-k^/c^  C  ’ 
1 


ct  + 


VI  -  fc2/c2  VI  - 


■  2:  . 


Transformation  (2.18)  introduces  the  new  coordinate  system  (r,  C,^).  The  origin  of  this  new  coordinate 
system  moves  with  the  speed  k  along  the  2:-axis  of  the  original  coordinate  system.  In  other  words,  at  every 
given  t  it  is  positioned  at  2;  =  kt  in  the  original  frame  of  reference.  In  implementing  transformation  (2.18), 
we  will  always  need  to  assume  that  k  <  c,  as  has  also  been  suggested  in  Section  2.2. 

The  key  property  of  the  Lorentz’  transform  (2.18)  is  that  it  does  not  change  the  form  of  the  wave 
equation  (2.1)  (and  consequently,  (2.15)  and  (2.17)),  see,  e.g.,  [25].  As  such,  let  us  introduce  an  arbitrary 
function  of  time  x  -  x(^).  x(t}  =  0  for  t  <  0,  so  that  it  also  be  smooth  on  the  entire  E.  Next,  we  define 
p2  _  y.2  _j_  ^2  ^  and  then 

i>(r,  C,  ^  ~  (2.19a) 


becomes  a  solution  to  the  wave  equation  in  the  new  coordinates  (r,  (,  0).  Solution  (2.19a)  is  driven  by  a  point 
(^-type  source,  which  is  located  at  the  origin  {t*  =  0,  C  =  0}  modulated  in  time  by  the  function  xW- 
;)^'(0)  =  0,  this  solution  also  satisfies  the  homogeneous  initial  conditions.  Consequently,  the  function 


p(r,z,t) 


(2.19b) 


obtained  by  substituting  (2.18)  into  (2.19a)  is  a  solution  of  equation  (2.15)  with  the  RHS  f{r,z^t)  =  x(^)  * 
6(r,z  —  kt).  In  other  words,  i{j{ryZjt)  of  (2.19b)  is  a  solution  to  the  wave  equation  excited  by  a  (J-source 
that  moves  straightforwardly  and  uniformly  and  is  modulated  in  time  by  a  given  smooth  function.  Solution 
(2,19b)  also  satisfies  homogeneous  initial  conditions  (2,2).  From  the  standpoint  of  physics,  solution  (2.19b) 
can  be  characterized  as  radiation  of  spherical  waves  by  a  moving  point  source. 

Solution  (2.19b)  is  obviously  singular.  To  use  it  for  testing  the  numerical  algorithm  we  need  to  remove 
the  singularity.  For  that,  let  us  first  introduce  the  actual  domain  S{t).  In  all  the  experiments  that  follow, 
the  domain  S{t)  is  a  sphere  of  diameter  d  with  its  center  at  the  origin  of  the  new  coordinate  system: 
{r  —  0,  z  =  kt].  As  such,  this  spherical  domain  moves  uniformly  along  the  z-axis,  which  obviously  helps  us 
keep  the  axial  symmetry  intact.  As  has  been  mentioned,  the  speed  of  this  motion  is  “subsonic,”  A:  <  c,  which 
conforms  to  one  of  the  key  requirements  for  building  the  lacunae-based  algorithm  previous!}'  put  forward  in 
Section  2.2. 
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Let  us  also  define  {z  —  kt)^  and  introduce  the  function  Q  =  Q(r),  r  >  0,  such  that  Q(0)  =  0, 

Q{r)  =  1  for  f  >  /cd/2,  where  k  <  1,  and  also  ^  ^  =  0  for  m  =  1, 2, . . .  till  at  least  m  =  4. 

Then,  it  is  easy  to  see  that  the  function  <^(r, =  'ip{r,z^t)  •  Q{r)  is  regular  (continuous  and  bounded) 
everywhere.  Moreover,  it  is  easy  to  verify  by  direct  differentiation  that  the  same  is  true  for  the  function 
f(r,z^t)  z,t),  where  □  denotes  the  wave  operator,  i.e.,  the  left-hand  side  of  equation  (2.15).  We 

will  use  /(r,  t)  defined  this  way  as  the  source  function  for  equation  (2.15).  Clearly,  /(r,  z,  t)  may,  generally 
speaking,  differ  from  zero  only  on  the  ball  of  a  smaller  diameter  Kd  concentric  with  S{t).  Everywhere  else, 
i.e.,  for  f  >  Kd/2,  f{r^z^t)  =  0. 

Obviously,  the  solution  of  problem  (2.15),  (2.2)  driven  by  this  /(r,  z,t)  —  Oip{r,z,t)  is  the  foregoing 

(p{r,  z,  t)  —  ip{r,  z,  t)  •  Q{r)  .  (2.20) 


This  function  satisfies  the  non-homogeneous  wave  equation  with  the  RHS  /(r,z,t)  on  a  smaller  ball  of 
diameter  Kd  concentric  with  S{t).  Every w'here  else  it  is  a  solution  to  the  homogeneous  wave  equation 
because  it  coincides  with  ip{r,Zjt)  of  (2.19b).  Consequently,  <^(r, z,t)  of  (2.20)  can  be  interpreted  as  the 
radiation  of  waves  by  a  compactly  supported  moving  source  /(r,  z,^).  Numerically,  we  will  be  reproducing 
solution  ip{r,z,t)  given  by  (2.20)  on  the  domain  S{t)  using  finite-difference  methods. 

W'e  employ  three  different  explicit  central-difference  schemes  in  our  simulations.  In  all  three  cases  we 
construct  a  uniform  rectangular  grid  on  the  plane  (r, z):  n  =  ^  =  0,1,...  ,iVr,  hr  —  RfNr,  and 

Zm  -  mhz,  m  =  0,  ±1, . . .  ,  d=iV^,  =  Z/2Nz.  The  discrete  time  levels  are  tn  =  nr,  n  =  0,1, ... .  For  the 
cell- centered  second-order  scheme,  we  keep  the  values  of  the  unknown  function  ip  at  the  grid  nodes  in  the  z 
direction  and  at  mid-points  in  the  r  direction: 


^1^1/2, m 


/+l/2,m 


,  1 


f— - 

\r  1+1/2  hr 


n+1 


m3/2.n 


•  v^r+1/2,. 


-n 


^1+1(2, m  ^f-l/2,m 


+ 


(2.21a) 


V^/Vl/2,m4-l  ^^1+1/2, m  + 

P 


/+l/2,m  • 


Equations  (2.21a)  hold  for  all  /  >  0.  As  in  this  case  we  do  not  have  the  unknown  function  defined  on  the 
axis  of  symmetry,  and  the  closest  values  that  correspond  to  /  =  0  are  half-grid-size  away: 


scheme  for  /  =  0  is  obtained  by  simply  assuming  that 


rt~*Pl-l/2,rr 


t/=0 


=  0,  which  can  be  interpreted  as 


a  second-order  approximation  of  the  natural  condition  ^  =0.  This  immediately  yields  for  /  —  0: 


r=0 


‘Pm.n.- ^^1/2, m  +  Vl72) 

(2.21b) 

2  f  1  1  ‘Ps/2,m  ~  ‘Pl/2,m  ‘^1/2, m+1  “  ‘^^l/2,m  +  'fl/2,m-l 

hr  hi 


TT^i' 

^1/2 


■)  —  fl/2,m  ■ 


For  the  node-centered  second-order  scheme,  tp  is  taken  at  the  actual  grid  nodes,  and  for  /  >  0  we  have: 


\n  hr 


- fZ - - hr 


(2.22a) 


'fi?,m+l-M,m+‘PT,m-l  \  _ 


hi 


)=/i%- 
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To  obtain  the  scheme  on  the  axis  of  symmetry  /  0  in  this  case,  we  need  to  approximate  equation  (2.17).  For 
the  ^  derivative  in  this  equation  we  can  first  formally  write  expression 

obviously  reduces  to  because  of  the  symmetry:  and  consequently,  w’e 

obtain: 


n— 1 


(4* 


I  ~  ‘Po.m  ‘Po,m+l  ^‘Po,m  +  'Po,m-l 

hi 


=  fln 


(2.22b) 


The  last  scheme  is  the  node-centered  fourth-order  scheme.  More  precisely,  it  approximates  spatial  derivatives 
with  the  accuracy  Oih^  -I-  h*)  and  temporal  derivative  with  the  accuracy  0{t^).  For  1  >  1  we  have: 


n+1 


\Srihr 


'^1+1/2 


hr 


-  n-1/2 


iflm  -  n-1, 

hr 


ii  J_ 

3  Ti  2hr 


V?r+2,m  W.m  " 

ihr  2hr 


(2.23a) 


•f 


-‘Pl,m+2  +  16¥>”m+l  - 


For  /  =  1  we  have  r/_i  =  ro  =  0  and  consequently: 


_2 


/4  1  1 

11  1 

\  3  ri  hr 

hr  hr  \ 

1 

ICO 

1 

2/1,. 

-‘Plm+2  +  16v^?.,n+l  -  30<„,  -h  -  ^lm-2 

I2hl 


d- 


(2.23b) 


Finally,  for  1  =  0  we  again  have  to  approximate  equation  (2.17).  Using  symmetry  like  in  the  previous  case, 
we  arrive  at: 


n-1 


2  I  o  -2<^?,m  +  32<™-30^ 


12/i2 


"¥^0,m+2  +  l^V^O,m+l  +  l^V^0,m-l  V^0,m-2  \  _  rn 

ml  ;  "  * 


(2.23c) 


For  all  three  schemes,  (2.21),  (2.22),  and  (2.23),  setting  the  discrete  boundary  conditions  (2.16)  on  the  outer 
boundary  of  the  auxiliary  domain  [0,R]  x  [-Z/2,Z/2]  is  straightforward.  An  extra  boundary  condition  is 
needed  for  the  fourth-order  approximation.  As  it  basically  does  not  matter  what  boundary  conditions  we 
use  on  the  outer  boundary  of  the  auxiliary  domain  (see  Section  2.2),  we  simply  set  “  0  addition 

to  m  ~  R-^gsirding  the  time  step  r,  all  three  schemes  are  explicit  and  as  such,  there  is  a  Courant-type 
stability  constraint. 

As  has  been  mentioned,  we  present  the  results  of  numerical  computations  that  follow  in  order  to  corrob¬ 
orate  the  theoretical  design  properties  of  the  lacunae-based  algorithm,  i.e.,  to  show  the  temporally  uniform 
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grid  convergence  on  long  time  intervals.  For  that  purpose,  we  conduct  a  grid  refinement  study,  i.e.,  approxi¬ 
mate  the  exact  solution  (2.20)  on  a  sequence  of  successively  more  fine  grids.  In  so  doing,  the  time  step  r  for 
the  two  second-order  schemes  (2.21)  and  (2.22)  is  always  reduced  with  the  same  rate  as  the  corresponding 
spatial  sizes  hr  and  /i.;  and  the  time  step  r  for  the  fourth-order  scheme  (2,23)  is  reduced  twice  as  fast 
(i.e.,  by  a  factor  of  four  every  time  hr  and  hz  are  reduced  by  a  factor  of  two)  so  that  to  demonstrate  the 
fourth-order  overall  convergence  in  the  end.  The  computations  in  each  case  were  run  till  the  dimensionless 
time  t  —  200*d/c,  i.e.,  for  200  times  the  time  interval  required  for  a  wave  to  cross  the  domain.  This  certainly 
qualifies  as  “long  term”  from  the  standpoint  of  any  conceivable  application. 


Lacunae-based  (2,2)  cell-centered  scheme  Lacunae-based  (2,2)  node-centered  scheme 


(a)  The  second-order  scheme  (2.21)  (b)  The  second-order  scheme  (2.22) 

Fig.  2.2.  Grid  convergence  study  for  the  long-term  lacunae-based  integration  of  the  wave  equation. 


Lacunae-based  (2,4)  node-centered  scheme 


Fig.  2.3.  Same  as  Fig.  2.2  for  the  fourth-order  scheme  (2.23). 


The  actual  parameters  that  we  have  used  for 
our  simulations  are  the  following:  R  =  Z  =  n, 
d  =  1.8,  c  =  1,  A;  =  0.2,  k  =  0.8.  The  spatial  grid 
is  composed  of  square  cells:  Nr  =  Nz  and  conse¬ 
quently,  hr  =  hz  =  h.  The  actual  grid  dimensions 
Nr  X  2Nz  are:  64  x  128,  128  x  256,  and  256  x  512. 
The  temporal  partition  size  2T,  see  (2.8),  is  found 
from  (2.6)  assuming  that  Tint  =  .  ^he 

overlap  parameter  a  =  1/2.  The  functions  0(t) 
and  Q{r)  on  the  intervals  of  their  variation  from 
0  to  1  are  built  as  polynomials  of  degree  9  (with 
only  odd  powers  included),  which  guarantees  four 
continuous  derivatives  in  transition  to  the  constant 
(either  0  or  1).  The  function  x(^)  is  defined  as  fol¬ 
lows:  x(^)  =  (1  +  jsint)  (1  -  where  P{t) 
is,  again,  a  polynomial  of  degree  9  that  decays 


15 


smoothly  from  1  to  0  on  the  interval  [0,1]  (four  continuous  derivatives).  Finally,  to  subtract  every  (pj 
from  the  overall  solution  at  the  proper  moment  t  =  {j  -  1  (7j)T  +  Tint,  we  first  march  equation  (2.10) 

from  t  =  (j  -  1  +  (yj)T  till  t  =  ( j  +  1  +  aj)T  and  then  use  Fourier  expansion  in  ^  and  expansion  with 
respect  to  the  corresponding  discrete  eigenfunctions  (calculated  numerically)  in  r  to  advance  it  further  till 
aj)T  +  Tint. 

In  Figure  2.2  we  show  error  profiles  (more  precisely,  natural  logarithm  of  the  relative  error  on  the  domain 
S{t)  in  the  maximum  norm  as  it  depends  on  the  dimensionless  time)  on  all  three  grids  for  both  second-order 
schemes  (2.21)  and  (2.22).  In  Figure  2.3,  similar  curves  are  shown  for  the  fourth-order  scheme  (2.23).  From 
these  figures  we  conclude  that  indeed  no  error  is  accumulated  in  the  course  of  computations  because  all  error 
profiles  are  flat  throughout  the  entire  200  •  d/c  time  interval.  Thus,  the  solution  does  not  deteriorate  as 
time  elapses.  Figure  2.2  also  shows  that  every  time  the  grid  is  refined  by  a  factor  of  two  the  error  drops  by 
approximately  a  factor  of  four,  which  indicates  the  second-order  convergence.  Similarly,  Figure  2.3  shows  that 
every  time  the  grid  is  refined  by  a  factor  of  two  the  error  drops  by  approximately  a  factor  of  sixteen,  which 
is  an  indication  of  the  fourth-order  convergence.  Consequently,  we  can  conclude  that  numerical  experiments 
fully  corroborate  the  theoretical  design  properties  of  the  algorithm. 

3.  Lacunae-Based  ABCs  for 
the  Wave  Equation.  The  lacunae- 
based  algorithm  of  Section  2  provides 
a  venue  for  constructing  the  ABCs  for 
a  class  of  problems  that  reduce  to  the 
homogeneous  wave  equation  in  the  far 
field.  We  schematically  depict  the  ge¬ 
ometric  setup  for  one  such  problem 
in  Figure  3.1,  assuming  for  simplic¬ 
ity  that  there  is  no  source  motion, 
A:  =  0,  and  the  computational  domain 
is  stationary.  We  emphasize  though 
that  this  is  not  a  limitation,  and  that 
the  actual  ABCs  will  be  constructed 
and  tested  for  the  general  case  of  a 
moving  computational  domain,  while 
the  law  of  motion  can  be  arbitrary, 
see  Section  2.1.  The  problem  to  be 

solved  on  the  bounded  interior  do- 
FIG.  3.1.  Schematic  geometric  setup  for  the  ABCs.  gee  Fig- 

ure  3.1,  may  involve  some  complex 
phenomena  whose  nature,  however,  is  not  essential  for  the  current  discussion.^  We  only  require  that  the 
overall  combined  formulation  of  the  problem  be  uniquely  solvable  and  well  posed  under  the  assumption  of 
radiation  of  waves  in  the  far  field  (from  S{t)  toward  infinity),  where  the  problem  is  assumed  to  be  governed 
by  the  homogeneous  w^ave  equation.  The  role  of  the  ABCs  (as  mentioned  in  Section  1)  is  to  provide  a  closure 
for  the  truncated  problem  solved  on  the  actual  computational  domain  S{t).  This  closure  has  to  ensure  that 

"^The  interior  domain  is,  of  course,  the  same  as  S{t)  of  Section  2;  for  the  stationary  case  we  obviously  have  S{t)  =  S{0). 
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the  corresponding  finite- domain  solution  recovered  with  the  help  of  the  ABCs  be  close  to  (ideally,  exactly 
the  same  as)  the  solution  of  the  original  non-truncated  problem  restricted  to  the  bounded  domain,  see  [2]. 

In  the  wave  propagation  framework  adopted  in  this  paper,  one  can  say  that  the  ABCs  have  to  replace 
the  entire  far  field,  i.e.,  everything  beyond  the  bounded  interior  domain  5(f),  so  that  the  resulting  artificial 
boundary  be  completely  transparent  for  all  the  outgoing  (i.e.,  radiated)  waves.  We  also  note  that  the 
incoming  waves,  provided  that  they  are  meaningful  for  a  particular  setup,  can,  in  fact,  be  taken  into  account 
through  the  boundary  conditions  as  well,  but  we  do  not  discuss  this  issue  here  for  the  reason  of  simplicity. 

3.1.  Preliminary  Considerations  in  the  Continuous  Framework.  Let  (pc  ^  (pc{Xjt)  he  3,  solution 
to  the  aforementioned  combined  problem.  In  the  far  field,  i.e.,  outside  5(f),  the  function  pc{x,  f)  satisfies  the 
homogeneous  wave  equation.  We  also  assume  for  simplicity  that  the  solution  (pc{x^t)  “smoothly  originates 
from  zero”  at  f  =  0,  (i.e.,  turns  into  zero  along  with  its  first  derivative)  in  much  the  same  way  as  the 
solution  <p{x,t)  of  (2.1),  (2.2)  does.  This  assumption,  in  fact,  will  present  no  limitation  when  constructing 
the  ABCs.  The  argument  is  the  same  as  the  one  that  allows  us  to  relax  the  assumption  of  homogeneity  of 
initial  conditions  when  building  the  original  lacunae-based  algorithm,  see  [23]. 

Let  us  now  introduce  a  special  multiplier  function  that  is  again  schematically  shown  in  Figure  3.1.  This 
function  //  =  /i(x,f)  is  defined  for  all  those  x  and  f,  for  which  the  solution  pc(x,t)  makes  sense.  We  first 
require  that  Vf  >  0,Vaj  ^  5(f)  :  /i(x,f)  =  1,  or  in  other  words,  that  the  multiplier  be  identically  equal  to 
one  everywhere  outside  the  computational  domain  5(f)  for  all  times.  We  also  require  that  the  multiplier 
be  identically  equal  to  zero,  fi(x,t)  =  0,  on  most  of  the  domain  5(f)  (again,  for  every  f)  except  next  to 
its  boundary  from  the  interior  side.  An  example  of  the  narrow  near-boundary  transition  region,  where 
the  multiplier  ji{x^t)  changes  its  value  from  zero  to  one,  is  shaded  on  Figure  3.1.  What  is  important,  we 
require  that  the  multiplier  /x(a?,f)  be  a  sufficiently  smooth  function  with  respect  to  both  x  and  f,  which 
essentially  means  that  the  transition  within  the  shaded  region  on  Figure  3.1  has  to  be  smooth.  Regarding 
the  time  dependency  of  f),  once  the  domain  5(f)  moves  according  to  a  prescribed  law,  the  construction 
of  the  multiplier  has  to  trace  that  motion.  If  the  computational  domain  is  stationary,  5(f)  =  5(0),  then  the 
multiplier  still  may,  but  does  not  have  to,  depend  on  time. 

Next,  we  apply  the  wave  operator  □  =  ^  —  c^A  of  (2.1)  to  the  function  /i(x,f)  •  pc{x,t),  which  is 
defined  everywhere,  i.e.,  both  inside  and  outside  5(f).^  We  will  obviously  have: 


□(m)  =  =  9{x,t) 


✓ 

=  0 

=  0 


Vf,  Vx  ^  5(f) 

in  the  transition  region 

“well  inside”  5(f) 


(3.1) 


The  function  g{x,t)  of  (3.1)  may  generally  speaking  differ  from  zero  only  in  the  foregoing  near-boundary 
transition  region;  it  is  zero  outside  5(f)  because  the  function  picpc  coincides  there  with  the  solution  pc 
the  homogeneous  wave  equation;  it  is  also  zero  inside  5(f)  because  ju  =  0  there.  On  Figure  3.1  the  non- zero 
portion  of  g{Xyt)  is  identified  as  the  right-hand  side  RHS. 

We  can  now  consider  the  problem  (2.1),  (2.2)  wnth  the  function  g{x^t)  of  (3.1)  substituted  instead  of  the 
generic  RHS  f(x,  t).  The  key  fact  that  we  will  need  for  constructing  the  ABCs,  and  that  follows  immediately 
from  the  unique  solvability  of  the  Cauchy  problem  for  the  wave  equation,  is  that  the  solution  to  this  problem 
will  coincide  with  /i(x,f)  •  pc{x^t)  everywhere.  What  will  be  of  particular  importance  to  us  is  that  as  such, 


®Note,  the  solution  may  not  be  defined  on  all  of  S{t)  if,  e.g.,  there  is  a  scatterer  inside.  As,  however,  —  0 

there,  w'e  can  consider  }x{x,t)  •  cpcix^t)  to  be  defined  everywhere. 
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this  solution  will  coincide  with  ipc{x^t)  outside  S{t)  for  all  times,  because  iJ>{x^t)  =  1  there.  In  other  words, 
we  have  replaced  all  of  the  interior  problem  on  S{t)  (no  matter  how  complex  it  may  be)  by  the  special 
near-boundary  source  function  g{x,t)  so  that  the  corresponding  far-field  portion  of  the  solution,  i.e.,  the 
wave-radiation  solution  outside  5(t),  remains  totally  unaffected.  Later  on,  see  Sections  3.2  and  3.3,  this 
reduction  interpreted  in  the  discrete  framework  will  be  used  for  setting  the  ABCs.  The  idea  is  to  use  the 
exterior  solution  obtained  in  an  alternative  way  through  integrating  the  near-boundary  sources  as  a  closure 
for  the  interior  problem  solved  on  the  finite  computational  domain. 


3.2.  The  Concept  of  Discrete  ABCs.  To  construct  the  ABCs  for  a  finite-difference  scheme  that 
approximates  the  problem  described  in  the  beginning  of  Section  3,  we  will  employ  the  considerations  similar 
to  those  of  Section  3.1,  but  on  the  discrete  level.  As  a  helpful  illustration,  we  will  first  consider  here  a 
one-dimensional  model  example,^  and  then,  in  the  next  Section  3.3,  show”  how  to  build  the  ABCs  for  the 
actual  multidimensional  wave  propagation  problems. 
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Assume  that  we  are  solving  a  one¬ 
dimensional  combined  problem  on  the  entire 
E.  The  computational  domain  S{t)  =  5  (i.e., 
the  “near  field”)  is  fixed,  it  is  the  half-line  x  < 
0  (more  precisely,  it  is  {{x,t)\x  <  0,  t  >  0}); 
its  complement  E\5  represents  the  “far  field,” 
which  is  to  be  truncated  and  replaced  with  the 
ABCs.  As  such,  the  ABCs  are  to  be  set  at  the 
interface  x  =  0.  In  accordance  with  the  pre¬ 
vious  discussion,  we  also  assume  that  the  far 
field  is  governed  by  the  one-dimensional  ho¬ 
mogeneous  wave  equation 


av 

dt^ 


(3.2a) 


which  is  approximated  by  the  standard 
second-order  central-difference  scheme 


Fig.  3.2.  Illustration  to  the  one- dimensional  example. 


C  -  _u. 


(3-2b) 


constructed  on  the  rectangular  grid  of  variables  x  and  t  with  sizes  h  and  r  ^  hfc  respectively,  using  the 
five- node  stencil  shown  in  Figure  3.2.  (Note,  all  the  schemes  used  for  simulations  in  Section  2.5,  see  (2.21), 
(2.22),  and  (2.23),  are  of  the  same  central-difference  explicit  three-level  type.  This,  however,  is  by  no  means 
a  limitation  —  the  ABCs  can  be  constructed  for  any  type  of  discretization.) 

To  create  the  discrete  near-boundary  sources  similar  to  those  of  Section  3.1,  and  eventually  set  the 
discrete  ABCs,  we  will  need  to  be  able  to  apply  inside  S  the  same  finite-difference  wave  operator  of  (3.2b) 
as  the  one  wre  are  using  on  E\iS'.  As  such,  we  formally  extend  the  exterior  discretization,  i.e.,  the  rectangular 
grid  with  h  x  r  ceils,  into  the  interior  domain  5,  as  shown  in  Figure  3.2.  We  re-emphasize,  however,  that 


® Generally  speaking,  one-dimensional  problems  do  not  have  lacunae  (except  in  special  cases);  as  such,  this  example  will 
only  demonstrate  the  formal  construction  of  the  ABCs  on  the  grid. 
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this  is  done  only  for  the  “artificial”  purpose  of  building  the  ABCs.  The  actual  governing  equation  in  the 
near  field,  i.e.,  on  5,  as  well  as  its  discrete  counterpart,  may  be  more  complex  than  the  wave  equations  (3.2) 
above.  In  fact,  neither  the  scheme  stencil  nor  the  grid  used  for  computations  in  the  near  field  have  to  be 
the  same  as  those  in  the  far  field  (although,  they,  of  course,  may).  We  only  require  that  the  exterior  scheme 
(3.2b)  wdth  the  stencil  shown  on  Figure  3.2  be  applicable  till  x  =  0.  And  all  we  need  to  know  from  the 
standpoint  of  setting  the  ABCs,  is  that  the  overall  combined  problem  (near  field  and  far  field)  is  uniquely 
solvable  and  well-posed. 

Let  us  now  consider  all  nodes  of  the  aforementioned  rectangular  grid  that  belong  to  S  (on  all  time  levels) , 
i.e.,  those  for  which  xj  ~  jh  <  0.  We  denote  this  set  of  grid  nodes  by  the  complementary  set  that 
consists  of  all  nodes  from  R\S,  i.e.,  those,  for  which  xj  >  0,  is  denoted  by  AT",  see  Figure  3.2.  If  we  formally 
apply  the  five-node  stencil  of  the  scheme  (3.2b)  to  each  node  from  A/’'^,  then  this  stencil  is  obviously  going 
to  sweep  one  more  vertical  row  of  nodes,  which  already  belongs  to  A/*”  (i.e.,  to  M\5),  and  which  is  denoted 
by  7“  on  Figure  3.2.  Reciprocally,  if  we  apply  the  stencil  to  every  node  in  A/*”,  it  will  also  sweep  the  nodes 
7+  that  are  already  in  A/*"^.  The  two-layer  grid  structure  7  =  7+  U  7“  will  be  called  the  grid  boundary;  it 
represents  on  the  discrete  level  the  continuous  interface  betw^een  S  and  ®\5,  which  is  the  vertical  line  x  =  0. 

Next,  we  assume  that  we  integrate  the  interior  problem  one  time  step  after  another,  and  that  we  already 
know’  the  discrete  solution  on  the  domain  5,  as  well  as  the  values  of  on  the  grid  boundary  7,  up  to 
a  certain  time  level  n  (in  particular,  n  may  be  equal  to  zero,  which  corresponds  to  the  initial  conditions).^ 
These  data  obviously  allow  us  to  advance  the  next  time  step  n  +  1  on  7“^  and  everywhere  inside  5;  in  so 
doing,  the  outermost  interior  location  7“*"  on  the  level  n  -h  1  is  computed  by  scheme  (3.2b)  using  the  stencil 
shown  in  Figure  3.2.  These  data,  however,  already  do  not  allow  us  to  calculate  the  discrete  solution 
at  7~  on  the  level  n  -f  1.  And  if  this  value  is  not  known,  then  we  cannot  advance  further  to  level 
n  +  2.  Therefore,  we  conclude  that  the  function  of  the  ABCs  in  the  discrete  framework  will  be  to  provide 
the  missing  boundary  values  of  the  solution  at  7“  on  all  time  levels,  one  after  another,  starting  from  n  =  1. 
This  indeed  constitutes  the  closure  of  the  discrete  system  solved  on  S. 

To  provide  the  foregoing  missing  boundary  value  for  a  given  n,  we  recall  that  even  so  we  do  not 
know  the  discrete  solution  on  level  n  -h  1  beyond  (i.e,  we  do  not  know  we  do  know  that  the 

solution  can  be  complemented  on  M"  to  a  solution  of  equation  (3.2b)  on  all  time  levels  till  n  -h  1.  For 
our  purposes,  we  will  only  need  the  existence  of  this  complement  rather  than  its  actual  representation.  Let 
us  now’  introduce  a  multiplier  function  p  similar  to  that  we  have  used  in  Section  3.1.  The  near-boundary 
interior  transition  region  for  this  multiplier  is  schematically  shown  by  the  shaded  area  on  Figure  3.2.  We 
apply  this  multiplier  to  the  combined  discrete  solution  levels  including  n  -f  1 

(obtaining  may  require  projecting  onto  the  grid  A/*'^).  In  so  doing,  nothing  changes  on  A/*“  U  7“^, 
because  //  =  1  for  x  >  0.  All  the  changes  due  to  multiplication  of  by  /x  will  obviously  be  introduced 
on  A^+\7'^  only.  Those  amount  to  a  smooth  passage  within  the  transition  region  (see  Figure  3.2)  from  the 
actual  unaltered  values  of  the  solution  on  7"*"  to  zero  “well  inside”  the  computational  domain. 

Next,  similarly  to  Section  3.1,  we  apply  the  discrete  wave  operator  of  (3.2b)  to  the  modified  solution 
As  is  defined  up  to  the  level  n  -1- 1,  the  result  will  be  defined  up  to  the  level  n.  Analogously 


■^We  use  the  subscript  “S'”  in  rather  than  “A/'+”  to  emphasize  that  the  actual  interior  discrete  solution  may  be  computed 
on  a  different  grid. 
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to  (3.1),  we  obtain  for  all  levels  till  n: 


=  gW 


=  0 
5^0 
=  0 


on  the  grid  A/*” 

on  7+  &  in  the  transition  region 
“well  inside” 


(3.3) 


Notice,  we  can  claim  that  the  result  in  (3.3)  is  zero  on  J\f~  only  because  nothing  has  been  modified  by  /x  on 
7"^  and  beyond.  As  such,  we  are  simply  using  the  fact  that  (pc^^  is  a  solution  to  the  homogeneous  equation 
outside  the  computational  domain,  and  we  do  not  need  to  know  the  explicit  form  of  this  solution. 

Finally,  we  solve  the  non-homogeneous  counterpart  to  equation  (3.2b)  driven  by  the  RHS  of  (3.3) 
everywhere  on  U  AT";  this  will  be  henceforth  referred  to  as  solving  the  auxiliary  problem.  According 
to  our  construction,  solving  the  auxiliary  problem  will  allow  us,  in  particular,  to  recover  the  value  of  , 
which  was  not  known  previously,  and  which  can  now  be  supplied  to  the  interior  scheme  as  the  missing 
boundary  value.  This  means  that  we  will  have  provided  the  ABCs  for  the  interior  problem,  because  in  so 
doing  we  complete  the  time  level  n  +  1  and  facilitate  advancing  the  next  level  n  +  2. 

The  are  several  important  comments  to  be  made  regarding  the  foregoing  ABCs’  algorithm.  At  a  first 
glance,  the  new  formulation  simply  does  not  change  much  from  the  standpoint  of  solving  the  original  infinite- 
domain  problem.  Indeed,  all  we  have  done  is  replaced  the  interior  problem  by  the  artificial  near-boundary 
sources  so  that  the  exterior  solution  remain  unaffected.  Then,  we  suggested  to  use  this  exterior  solution 
to  close  the  interior  discretization.  However,  obtaining  this  solution,  i.e.,  solving  the  auxiliary  problem, 
basically  brings  along  the  exact  same  set  of  complications  that  we  have  been  trying  to  avoid  by  introducing 
the  ABCs.  Indeed,  the  domain  of  the  auxiliary  problem  A/*"*"  U  A/*”  is  still  unbounded  and  as  such,  special 
treatment  will  be  required  for  its  numerical  solution. 

There  is,  however,  a  fundamental  difference.  The  new  auxiliary  problem  is  linear  throughout  the  entire 
space,  and  it  is  driven  by  known  sources  that  are  compactly  supported  inside  the  computational  domain  5. 
Consequently,  the  lacunae-based  algorithm  of  Section  2  appears  to  be  a  most  natural  tool  to  solve  it.^  Em¬ 
ploying  the  lacunae-based  algorithm  immediately  implies  that  the  domain  of  the  auxiliary  problem  becomes 
bounded.  Moreover,  the  “sufficiently  retarded”  sources  do  not  contribute  to  its  solution  (see  Section  2),  i.e., 
only  a  limited  extent  of  temporal  pre-history  of  the  solution  will  be  needed  to  sustain  the  continuous  time 
marching  no  matter  how  far  in  time  we  would  like  to  advance  the  solution.  In  other  words,  the  missing 
boundary  value  for  the  interior  discretization  can  be  obtained  using  only  finite  computer  resources  in 
terms  of  both  memory  and  number  of  arithmetic  operations.  Furthermore,  these  resources  (say,  per  time 
level)  will  not  increase  no  matter  for  how  long  we  may  need  to  run  the  computation,  i.e.,  how  large  n  may 
become.  In  this  sense,  the  proposed  ABCs  become  “true  ABCs,”  i.e.,  the  procedure  that  guarantees  the 
appropriate  closure  of  the  truncated  problem  with  only  finite  non-growing  amount  of  computer  resources 
required.  In  addition  to  that,  we  are  guaranteed  that  the  ABCs  as  a  part  of  the  overall  algorithm  will  not 
contribute  toward  the  buildup  of  numerical  error  during  long  runs. 

The  proposed  ABCs  can  obviously  be  implemented  via  alternating  interior /exterior  steps.  Namely,  we 
advance  one  time  step  in  the  interior  (including  7+)  assuming  that  all  the  data  that  we  need  from  the 
previous  time  levels  are  available.  The  resulting  newly  calculated  time  level  will  be  the  only  one  to  which  the 
multiplier  has  not  been  applied  yet.  We  multiply  it  by  /x,  and  then  apply  the  direct  operator  thus  obtaining 


®We  reiterate  that  this  algorithm  cannot  be  applied  in  the  case  of  one  space  dimension,  but  our  ultimate  goal  is  three- 
dimensional  problems,  see  Section  3.3,  and  the  considerations  of  the  current  section  are  for  illustrative  purposes  only. 


20 


the  right-hand  side  see  (3.3),  on  one  more  time  level  as  well.  Finally,  we  perform  one  step  of  the 
lacunae-based  integration  of  the  auxiliary  problem  driven  by  and  obtain  the  missing  boundary  value  for 
the  interior  problem.  Then,  the  procedure  cyclically  repeats  itself.  Summarizing,  we  can  say  that  having 
advanced  the  interior  solution,  we  then  generate  a  new  contribution  to  the  RHS  of  the  auxiliary  problem 
and  subsequently  advance  its  solution,  which,  in  turn,  allows  us  to  calculate  the  next  interior  step. 

An  important  observation,  which  is  easy  to  make,  is  that  the  missing  boundary  value  that  the 
ABCs  provide  does  not,  of  course,  depend  on  the  actual  shape  of  the  multiplier  ju  in  the  transition  region. 
Indeed,  ji  is  defined  so  that  it  does  not  alter  the  solution  on  the  grid  boundary  7.  Consequently,  when  we 
first  apply  the  direct  operator  to  see  (3.3),  and  then  integrate  the  non-homogeneous  wave  equation 

driven  by  g^^\  the  solution  on  7  will  remain  unchanged  no  matter  what  changes  have  been  introduced  by  /i 
in  the  interior.  As  such,  the  value  will  only  depend  on  the  values  of  the  solution  on  7  on  all  previous 
time  levels,  as  well  as  on  Moreover,  since  all  the  operations  that  we  perform  when  constructing  the 

ABCs  are  linear,  we  can  symbolically  write  the  resulting  boundary  condition  as  a  linear  form: 


n+l 

7“ 


(3.4) 


Technically,  the  dependency  of  on  the  previous  time  levels,  see  (3.4),  involves  all  of  the  latter,  from 
all  the  way  back  till  n  =  0.  However,  the  use  of  the  lacunae  in  three  space  dimensions  will  allow  us 
to  truncate  (3.4)  and  leave  only  several  levels  that  immediately  precede  n  -I- 1;  the  number  of  the  levels 
involved  will  be  fixed  and  will  not  increase  with  the  increase  of  n.  As  such,  temporal  non-locality  of  the 
ABCs  will  be  limited,  and  this  will  not  be  a  consequence  of  any  approximation,  but  rather  an  implication  of 
the  fundamental  properties  of  the  problem.  We  also  note  that  the  representation  of  the  ABCs  in  the  form 
(3.4)  serves  primarily  the  reason  of  convenience  and  compactness  in  notations.  In  fact,  the  coefficients  of 
the  linear  form  I  never  need  to  be  known  explicitly  except,  possibly,  when  multiple  interior  problems  are 
solved  with  the  same  exterior  model  (means  same  grid,  same  geometry,  same  scheme).  In  this  case  it  may 
be  beneficial  to  calculate  the  form  I  once  and  ahead  of  time,  compared  to  the  straightforward  calculation  of 
many  times  according  to  the  procedure  outlined  above. 

It  is  also  important  to  mention  that  boundary  condition  (3.4)  can  be  obtained  in  the  framework  of  a 
general  unsteady  ABCs’  methodology  proposed  by  Ryaben’kii  in  [24]  (see  also  older  work  [26])  for  a  variety 
of  problems,  including  multidimensional  cases,  domains  of  varying  shape,  and  different  types  of  schemes 
explicit  as  well  as  implicit.  Work  [24,26]  describes  the  theoretical  construction  of  the  ABCs  per  se,  and  does 
not  address  any  issues  related  to  the  actual  computations  (for  example,  using  lacunae-based  integration,  as 
proposed  in  the  current  paper).  The  methodology  of  [24,26]  relies  on  the  concepts  of  generalized  potentials 
and  boundary  projection  operators  of  Calderon’s  type  obtained  and  implemented  in  the  discrete  framework 
by  means  of  the  difference  potentials  method,  see  [27-29].  In  this  perspective,  the  ABCs  of  [24,26],  and 
boundary  condition  (3.4)  in  particular,  can  be  interpreted  as  discrete  counterparts  to  Calderon’s  boundary 
equations  with  projections  in  the  unsteady  case.  In  the  following  Section  3.3  we  will  describe  a  direct 
approach  to  obtaining  multidimensional  ABCs  on  moving  boundaries,  with  no  explicit  use  of  the  apparatus 
of  Calderon’s  projections,  and  will  also  show  how  to  apply  the  lacunae-based  algorithm  to  perform  the 
computations  needed  for  these  boundary  conditions. 

To  conclude  this  section,  we  emphasize  that  even  so  obtained  according  to  (3.4)  formally  does  not 
depend  on  the  shape  of  the  multiplier  pi  inside  5,  we  still  need  to  have  this  multiplier  smooth.  In  other 
words,  we  could  not  have  used,  e.g.,  a  step  function,  instead  of  pi.  The  reason  is  that  non-smoothness  will 
ruin  lacunae  in  the  discrete  solution  (see  Section  2  and  [23]  for  more  detail)  and  consequently,  we  will  no 
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longer  be  able  to  use  the  lacunae-based  integration  for  solving  the  auxiliary  problem  and  as  such,  setting 
the  ABCs. 


3.3.  General  Construction  of  Discrete  ABCs.  Similarly  to  the  setup  of  Section  2,  we  now  consider 
a  domain  S{t)  €  that  has  finite  diameter  d  for  all  times,  and  other  than  that  may  travel  in  space  according 
to  a  prescribed  law  with  the  only  limitation  that  its  maximum  speed  be  subsonic:  k  <  c.  S{t)  will  be  the 
computational  domain,  or  near  field.  In  the  far  field,  i.e.,  outside  S{t),  we  assume  that  our  model  is  governed 
by  the  homogeneous  wave  equation: 


dt‘^  \dxl  dx\  dxl) 


=  0, 


t  >  0  . 


(3.5) 


As  we  have  discussed,  inside  S{t)  the  solution  ip  —  (p{x,t)  may  be  governed  by  a  more  complex  equa¬ 
tion/system,  but  all  we  need  to  assume  is  that  the  overall  problem  be  uniquely  solvable  and  well  posed  under 
the  condition  of  waves’  radiation  toward  infinity.  For  simplicity,  we  also  assume  homogeneity  of  the  initial 
data  everywhere,  which  is,  however,  not  a  limitation  (see  [23]), 

Let  us  now  introduce  the  discretization  grid  for  equation  (3.5).  In  principle,  we  need  this  grid  only  in  the 
far  field,  i.e.,  outside  5(t),  because  the  interior  problem  may  be  discretized  in  a  different  way,  as  indicated 
before.  As,  however,  we  have  also  seen,  to  obtain  the  ABCs  we  need  to  set  up  the  auxiliary  problem  for  the 
non-homogeneous  counterpart  of  equation  (3.5)  driven  by  the  special  near-boundary  sources.  The  auxiliary 
problem  is  to  be  formulated  and  solved  on  the  entire  space.  As  such,  we  introduce  the  grid  for  the  linear  wave 
equation  on  the  entire  x  [0,  oo)  as  well.  We  denote  by  Af  the  collection  of  all  grid  nodes  in  E^  x  [0,  oo), 
on  which  w'e  evaluate  the  solution  p.  Since  (3.5)  is  an  evolution  equation,  it  is  convenient  to  consider  Af  as 
a  composition  of  spatially  aligned  grid  hyper-planes:  J\f  =  AfoUAfiU...UAfn  —  Each  A/k  is  a  spatial  grid 
on  E^ ,  and  we  emphasize  that  they  may,  but  do  not  have  to,  be  the  same  on  different  levels  n. 

Let  the  individual  nodes  of  the  grid  Af  be  denoted  by  n.  Equation  (3.5)  is  approximated  by  a  finite- 
difference  scheme,  which  we  assume,  of  course,  to  be  consistent  and  stable: 


^  amnPn  =  0.  (3-^) 

neJ^m 

In  (3.6),  Afm  denotes  the  stencil  attributed  to  the  node  m,  and  Umn  are  the  corresponding  coefficients.  When 
we  say  that  the  stencil  is  attributed  to  a  particular  node,  we  mean  that  the  residuals  of  the  discrete  equation 
are  evaluated  at  this  particular  grid  location.  In  regard  to  this,  we  note  that  the  residuals  of  the  discretized 
equation  (3.5)  may,  but  do  not  have  to,  be  evaluated  on  the  same  grid  Af.  To  preserve  the  generality  of 
the  discussion,  w’^e  assume  that  there  is  another,  different,  grid  in  E^  x  [0,  oo),  on  which  we  keep  the 
residuals,  as  well  as  the  right-hand  sides,  if  any,  of  the  discrete  wave  equation.  The  subscript  m  in  equation 
(3.6)  basically  refers  to  this  grid:  m  e  M.  In  the  one-dimensional  example  of  Section  3.2,  both  grids  Af  and 
Ai  were  simply  the  same,  and  we  did  not  have  to  distinguish  between  the  two.  To  give  an  example  of  the 
opposite  type,  we  mention  the  Yee  scheme,  see  [30],  which  is  one  of  the  primary  tools  for  discretizing  the 
Maxw^ell’s  equations,®  and  which  involves  staggering  in  both  space  and  time. 

Next,  w’e  introduce  two  subsets  of  nodes  of  the  grid  Af.  Let  the  level  A/*n  correspond  to  the  actual  moment 
of  time  For  every  n,  we  define  Af^  as  the  set  of  all  those  and  only  nodes  on  this  level  that  belong  to  the 
domain  S{tn),  and  Af^  as  the  complement  of  A/*^  to  the  entire  Afn’  Af^  =  Afn\Afn-  In  other  words,  Afn 

^The  simplest  version  of  the  Maxwell’s  equations  describes  the  propagation  of  electromagnetic  waves  in  vacuum,  which  is 
a  w^ave  model  similar  in  many  respects  to  that  given  by  (3.5). 


22 


contains  all  those  and  only  those  nodes  of  jVn  that  belong  to  R^\S{tn).  Subsequently,  we  define  the  set 
as  the  composition  of  all  A/'+  for  all  levels,  and  the  set  J\f~  as  the  composition  of  all  Af~  for  all  levels: 

U  AC,  U  AC: 

n  n 

Clearly,  A^"  =  A/*\A/*'^. 

In  our  definition  of  the  scheme,  see  (3.6),  we  have  identified  the  stencil  AC  with  the  grid  location  m, 
at  which  the  residual  is  evaluated.  From  here  on,  we  will  assume  for  simplicity  that  the  scheme  (3.6)  is 
explicit.  In  this  case,  there  is  only  one  non-zero  coefficient  amn  on  the  upper  time  level  of  the  stencil  AC- 
We  will  denote  the  corresponding  grid  node  by  n,  and  when  it  may  not  cause  confusion,  we  will  refer  to  the 
same  stencil  as  either  AC  or  It  will  also  be  convenient  to  introduce  the  four-dimensional  (space-time) 
vector  b  =  h  —  m.  The  meaning  of  this  vector  is  that  is  defines  the  relative  position  of  the  node  n,  at  which 
the  upper-level  coefficient  is  non-zero:  a^n  ^  0?  with  respect  to  the  “center”  of  the  stencil  m.  This  vector 
is  obviously  constant,  it  depends  only  on  the  local  structure  of  the  stencil,  and  does  not  depend  on  where 
exactly  on  the  grid  this  stencil  is  applied  at  every  given  moment.  In  the  one-dimensional  example  considered 
previously,  we  would  have  b  =  (r,  0) . 

Let  us  now  apply  the  stencil  Afn  to 
every  node  n  G  AT^;  in  so  doing,  the 
stencil  obviously  sweeps  the  entire  grid 
A/^+,  as  well  as  a  portion  of  the  grid  A/*“ 
next  to  the  interface,  we  will  denote  this 
portion  by  7": 


On  Figure  3.3,  we  present  a  one¬ 
dimensional  illustration^®  similar  to 
that  in  Figure  3.2  for  the  case  of  a  uni¬ 
form  motion  of  the  computational  do¬ 
main,  when  the  space-time  trajectory  of 
the  boundary  is  a  straight  line  (the  set 
7“  is  denoted  by  small  circles).  From 

the  standpoint  of  implementation,  the 
Fig.  3.3.  One-dimensional  illustration  to  the  case  of  a  moving  domain.  _ 

values  of  the  solution  at  the  nodes  7 

are  exactly  those  that  need  to  be  provided  by  the  ABCs  from  the  exterior  side  so  that  to  be  able  to  calculate 
the  solution  at  every  interior  node  n  G  A/’“*'  using  the  scheme  (3.6).  Reciprocally,  the  stencil  Mh  applied  to 
every  node  h  G  sweeps  additional  nodes  7“^  C  A/''*',  see  “bullets”  on  Figure  3.3: 

The  set  7+  complements  7"  to  the  complete  grid  structure  known  as  the  grid  boundary  7  (see  [28]): 

7  =  7+U7“- 

note  again  that  everywhere  in  this  section  the  one-dimensional  examples  are  for  illustration  purposes  only,  the  actual 
algorithm  is  three-dimensional. 
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The  grid  boundary  is  a  multi-layer  fringe  of  nodes  (two-layer  fringe  in  the  particular  case  of  a  second-order 
scheme  with  the  stencil  depicted  in  Figure  3.3)  that  is  located  near  the  continuous  boundary  and  straddles 
it  in  some  sense  (cf.  to  Section  3.2). 

To  proceed  with  the  construction  of  the  ABCs,  we  will  need  to  assume  hereafter  that  the  region  of 
linearity,  i.e.,  the  area  where  the  solution  of  the  overall  combined  problem  is  governed  by  the  linear  homoge¬ 
neous  w’ave  equation,  extends  “a  little  bit”  to  the  interior  of  the  computational  domain  S{t)  as  well.  More 
precisely,  this  region  will  be  assumed  to  extend  inward  at  least  as  far  as  the  entire  grid  boundary  7.  This 
obviously  presents  no  limitation  from  any  standpoint.  The  multiplier  /i  =  M(a;,t)  in  this  case  is  required  to 
be  identically  equal  to  one,  =  1,  not  only  outside  S{t),  but  also  inside  —  again,  to  the  extent  of  7”^. 

As  such,  the  transition  region  for  the  multiplier,  which  is  schematically  shown  by  the  darker  gray  shading 
on  Figure  3.3,  is  shifted  away  from  the  boundary  of  S{t). 

To  actually  build  the  ABCs,  we  will  perform  the  same  procedure  as  outlined  in  Section  3.2.  What  we 
actually  need  in  the  discrete  framework  is  to  obtain  the  missing  exterior  boundary  values  of  the  solution 
on  every  time  level  n  -f  1,  or  in  other  words,  to  complete  this  time  level,  so  that  to  be  able  to  advance 
the  next  time  step.  To  do  that,  we  take  the  solution  already  computed  inside  S(t)  up  to  the  level  n  -h  1, 
multiply  it  by  /i  and  then  apply  the  discrete  operator  everywhere.  In  doing  so  we  assume,  as  mentioned 
before,  that  starting  from  7  outwmd,  the  solution  satisfies  the  discrete  homogeneous  wave  equation.  In  the 
general  case  that  we  are  looking  at  now,  an  application  of  the  operator  brings  us  from  the  grid  AT  to 
the  grid  A4 .  The  construction  of  the  grid  boundary  7  and  multiplier  fj,  guarantees  that  on  all  time  levels  up 
to  n  the  near-boundary  artificial  sources  will  satisfy  (cf.  to  (3.3)): 

=  0  for  such  m  £  M  that  m  +  6  e  AT" 

^  0  for  such  near-boundary  m  £  M  that  m  +  b  £  J\f~^  (3.7) 

—  0  one  the  grid  M  “well  inside”  S  (t) 

Formula  (3.7)  suggests  that  in  addition  to  the  actual  interface  between  the  domains,  i.e.,  the  boundary  of 
S(t),  it  will  also  be  convenient  to  consider  another  space-time  trajectory  obtained  from  this  original  interface 
by  the  constant  displacement  —6,  it  is  shown  by  the  dash-dotted  line  in  Figure  3.3.  The  right-hand  side  g^^^ 
will  be  zero  on  A4  everywhere  outside  this  new  displaced  boundary,  and  will  differ  from  zero  right  next  to  it 
on  the  interior  side.  On  Figure  3.3,  we  schematically  show  by  the  lighter  gray  shading  the  region  where  we 
still  have  ju  =  1  but  the  RHS  g^^^  may  already  differ  from  zero.  Note,  in  the  example  of  Section  3.2  be  did 
not  need  to  consider  the  displaced  interface  because  the  domain  S(t)  was  stationary  and  the  displacement 
b  was  parallel  to  the  time  axis:  b  =  (r,  0).  Another  important  case  when  the  two  boundaries  would 
coincide  isn  —  m  —  b  =  0.  However,  we  cannot  generally  assume  that  for  time-dependent  problems. 
On  the  other  hand,  we  mention  that  the  grid  boundaries  originally  introduced  in  [27,  28]  and  previous 
publications  by  Ryaben’kii  for  the  solution  of  steady-state  problems  using  difference  potentials  method,  have 
been  constructed  so  that  just  the  center  m  of  the  stencil  Mm  (where  the  residuals  are  evaluated)  would  sweep 
a  given  grid  subdomain  and  as  such,  generate  the  aforementioned  fringe  of  nodes  7  next  to  the  continuous 
boundary. 

Let  us  now  make  a  few  remarks  of  explanatory  nature  regarding  the  structure  of  the  grid  boundary  7. 
It  obviously  depends  only  on  the  type  of  the  stencil  Mm  and  geometry  of  the  actual  continuous  boundary 
that  it  straddles.  From  the  definition  of  7  is  is  easy  to  see  that  once  we  have  a  solution  to  the  homogeneous 
equation  on  7  and  every wdiere  in  the  exterior,  and  operate  by  on  this  solution,  then  we  can  guarantee 
without  actual  calculation  that  the  result  will  be  zero  for  all  m  €  A4:  m  +  b  £  A/*”.  However,  we  cannot 
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“touch”  even  one  single  node  from  7  so  that  not  to  lose  this  property.  If,  for  example,  we  allow  for  an 
alteration  (via  p  1)  of  a  node  from  7+  (see  Figure  3.3),  it  will  necessarily  affect  the  exterior  RHS.  The 
latter  may,  generally  speaking,  become  non-zero  at  some  node(s)  m  e  M:  m  +  5  €  A/*“,  and  we  will  no 
longer  be  able  to  actually  calculate  it  because  we  do  not  know  anything  about  the  exterior  solution  beyond  7 
except  that  it  satisfies  the  homogeneous  equation.  We  see,  therefore,  that  it  would  ruin  the  entire  derivation. 
On  the  other  hand,  the  construction  of  the  grid  boundary  7  is  consistent  in  the  sense  that  to  calculate  the 
actual  non-zero  near  boundary  sources  for  those  m  G  for  which  m  -f  6  G  it  is  sufficient  to  know 
the  solution  only  on  7  and  further  inward,  nothing  outside  7  needs  to  be  known.  As  for  the  values  that 
are  still  needed,  those  are  provided  by  the  ABCs’  algorithm  on  every  time  level  and  as  such,  are  available 
on  all  subsequent  levels  for  calculating  the  source  terms  of  the  auxiliary  problem. 

Having  outlined  the  construction  of  the  grid  boundary  7  and  near-boundary  sources  in  the  general 
case,  we  build  the  actual  ABCs’  algorithm  in  much  the  same  way  as  described  previously.  We  perform  the 
alternating  interior/exterior  steps:  First  advance  one  step  in  the  interior,  then  apply  fx  and  calculate  one 
more  level  of  the  sources  g^^^ ,  and  finally  make  one  step  of  the  lacunae-based  integration  of  the  auxiliary 
problem  driven  by  these  sources  (see  Section  2),  thereby  providing  the  missing  data  for  advancing  the  next 
interior  step.  Then,  the  procedure  cyclically  repeats  itself.  To  solve  the  auxiliary  problem  in  this  general  case, 
we  will  obviously  need  a  full-fledged  version  of  the  lacunae-based  algorithm  (see  Section  2)  that  accounts  for 
the  motion  of  the  sources  and  employs  periodic  boundary  conditions  in  space  and  continuous  time-marching 
with  cyclic  subtraction  of  the  retarded  contributions.  As  has  been  shown,  implementation  of  the  lacunae- 
based  integration  technique  guarantees  that  the  domain  of  the  auxiliary  problem  will  be  bounded,  and  the 
computer  resources  needed  for  the  ABCs  will  be  finite  and  will  not  grow  with  time.  Other  properties  of 
the  ABCs  outlined  in  Section  3.2,  namely,  independency  on  the  shape  of  the  multiplier  /i,  and  possibility  to 
express  the  boundary  values  on  the  current  time  level  as  a  linear  function  of  the  values  on  the  previous  levels, 
see  (3.4),  hold  in  the  general  framework  of  this  section  as  well.  Because  of  the  lacunae,  the  aforementioned 
linear  form  will  depend  only  on  the  finite  non-increasing  number  of  the  preceding  time  levels  n,  essentially 
those  included  in  the  summation  (2.11)  once  this  formula  is  discretized  on  the  grid  Af .  This  means  that  the 
temporal  nonlocality  of  the  ABCs  will  be  limited.  As  for  the  multiplier  /i,  it  has  to  be  chosen  sufficiently 
smooth  so  that  to  maintain  good  quality  of  the  lacunae  in  the  discrete  solution,  see  Section  4. 

4.  Numerical  Experiments  with  the  ABCs. 

4.1.  The  Wave  Equation  with  Known  Exact  Solution.  The  first  case  that  we  analyze  in  the 
framework  of  the  ABCs  is  actually  the  exact  same  problem  that  we  solved  in  Section  2.5.  It  is  the  wave 
equation  driven  by  a  compactly  supported  oscillatory  source  in  straightforward  uniform  motion.  The  exact 
solution  of  the  problem  was  obtained  using  the  Lorentz  transform.  The  key  difference  between  the  current 
approach  and  that  of  Section  2.5  is  that  previously  we  have  applied  the  lacunae-based  algorithm  directly  to 
the  original  problem.  Here,  we  rather  decompose  the  problem  into  the  near  field  S{t)  and  far  field  M^\5(t), 
even  so  both  are  governed  by  the  same  wave  equation.  The  integration  in  the  near  field  is  then  performed 
by  the  conventional  time  marching.  The  exterior  closure  needed  to  sustain  this  time  marching  is  provided 
by  the  discrete  ABCs  on  the  boundary  of  S{t).  The  ABCs  are  constructed  on  the  basis  of  the  procedure 
outlined  above  —  through  the  lacunae-based  integration  of  the  artificial  near-boundary  sources. 

Similarly  to  Section  2.5,  we  have  implemented  three  different  schemes  (2.21),  (2.22),  and  (2.23),  and 
every  time  integrated  the  problem  till  the  dimensionless  time  has  reached  t  =  200  •  d/c.  The  multiplier 
/X  =  pi{r,z,t)  was  constructed  so  that  to  have  four  continuous  derivatives  with  respect  to  all  its  arguments. 
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Smooth  transition  from  0  to  1  was  obtained  with  the  help  of  algebraic  polynomials  of  degree  9  similar  to 
those  used  in  Section  2.5.  The  extent  of  the  transition  region  has  varied  slightly  between  different  cases  with 
no  noticeable  effect  on  the  quality  of  the  solution.  For  all  computational  variants  that  we  have  considered  it 
was  within  the  range  of  several  grid  cells  (typically,  eight  to  ten,  see  Section  4.4  for  further  detail). 

Lacunae-based  ABCs  for  (2,2)  celhcentered  scheme  Lacunae-based  ABCs  for  (2,2)  node-centered  scheme 
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(a)  The  second-order  scheme  (2.21)  (b)  The  second-order  scheme  (2.22) 

Fig.  4.1.  Grid  convergence  study  with  the  lacunae-based  ABCs  for  the  wave  equation. 


Lacunae-based  ABCs  for  (2,4)  node-centered  scheme 


Dimensionless  time 


Fig.  4.2.  Same  as  Fig.  4-^  for  the  fourth-order  scheme  (2.23). 


In  Figures  4.1  and  4.2  we  present  the  results 
of  the  grid  convergence  study  for  the  wave  equa¬ 
tion  integrated  with  the  lacunae-based  ABCs  over 
a  long  time  interval.  The  errors  are  evaluated  on 
the  interior  domain  S{t)  in  the  maximum  norm. 
Computational  setup  completely  follows  that  of 
Section  2.5,  where  we  applied  the  lacunae-based 
algorithm  to  the  original  problem  directly,  see  Fig¬ 
ures  2.2  and  2.3.  Namely,  each  scheme  was  im¬ 
plemented  on  a  sequence  of  three  grids:  64  x  128, 
128  X  256,  and  256  x  512;  all  geometric  parameters, 
grid  sizes,  parameters  that  control  the  partition 
(2.7),  as  well  as  the  actual  exact  solution  (2.19), 
(2.20),  against  which  we  compare  our  numerical 
results,  were  taken  exactly  the  same  as  before.  We 
only  note  that  in  this  case  partition  (2.7)  applies  to 


the  artificial  near-boundary  sources  needed  for  constructing  of  the  ABCs,  rather  than  the  original  right-hand 


side  /  that  drives  equation  (2.15). 

An  obvious  observation  w^hich  is  easy  to  make  is  that  Figures  4.1(a),  4.1(b),  and  4.2,  look  practically 
indistinguishable  from  Figures  2.2(a),  2.2(b),  and  2.3,  respectively.  In  other  words,  the  actual  levels  of  the 
error  on  corresponding  grids  are  essentially  the  same.  As  such,  we  conclude  that  in  this  most  simple  case  the 
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introduction  of  the  ABCs  makes  the  outer  boundary  of  the  computational  domain  completely  transparent 
for  all  the  outgoing  waves.  This  is  equivalent  to  saying  that  the  external  boundary  generates  no  reflection 
or  alternatively,  that  any  imperfections  associated  with  the  treatment  of  the  outer  boundary  can  always  be 
kept  on  or  below  the  level  of  the  truncation  error  pertinent  to  the  interior  discretization.  In  this  sense,  the 
discrete  lacunae-based  ABCs  that  we  have  constructed  can  be  regarded  as  an  ideal  closure  of  the  interior 
finite-difference  scheme.  Experimentally,  this  is  corroborated  by  the  fact  of  non- deteriorating  convergence  of 
the  scheme  with  the  theoretically  prescribed  rate  to  the  specially  constructed  exact  solution  of  wave-radiation 
type  on  the  computational  domain  S{t), 

4.2.  Nonuniform  Motion  of  the  Source.  The  exact  solution  (2.19)  of  the  wave  equation  driven  by 
a  point  source  in  straightforward  uniform  motion  was  obtained  in  Section  2.5  using  Lorentz’  transform.  The 
same  solution  could,  of  course,  have  been  obtained  by  directly  applying  the  Kirchhoff  integral  (2.3).  The 
integration  would  require  calculating  explicitly  the  location,  at  which  the  trajectory  of  the  source  intersects 
the  characteristic  cone,  and  would  lead  to  the  same  analytic  result  through  a  somewhat  more  complex 
derivation.  In  this  section,  we  consider  a  somewhat  more  complex  case  of  straightforward  nonuniform  (i.e., 
accelerated)  motion  of  the  source. The  Lorentz  transform  will  obviously  not  apply  in  this  case,  but  the 
Kirchhoff  integral  can  still  be  used  for  obtaining  the  exact  solution.  However,  for  the  general  accelerated 
motion  it  may  be  impossible  to  analytically  find  the  intersection  of  the  characteristic  cone  with  the  source 
trajectory.  Basically,  it  will  require  numerical  computation,  thus  making  the  resulting  exact  solution  only 
“semi- analytic.”  As  such,  to  analyze  the  case  of  the  accelerated  motion,  we  have  rather  chosen  a  full- 
fledged  numerical  approach.  First  we  calculate  a  fine-grid  reference  solution  using  the  original  lacunae-based 
algorithm  of  Section  2,  and  then  compare  against  it  the  solutions  obtained  with  the  ABCs  on  coarser  grids. 

Taking  the  value  of  the  speed  of  sound  to  be  equal  to  one,  c  —  1,  we  consider  the  following  law  of 
accelerated  motion  for  the  center  of  the  domain  S{t): 

r  =  0,  2;  =  zo{t)  =  kt-{-  k{cost  -  1), 

where  k  —  0.1.  Keeping  the  values  of  all  geometric  parameters  (sizes  of  the  domains,  etc.)  the  same  as 
before,  see  Section  2.5,  w^e  introduce  the  following  excitation  for  the  wave  equation  (2.15): 

Hr. =  COS  (II)  ^  P  (I)  .  (l  +  i  si„(V2.))  ■  P  (1  -  I;)  .  (4.1) 

where  d  =  1.8  ac  =  0.8,  =  p-\-{z-  ^o(^))^  and  P(-)  is  the  polynomial  of  9th  degree  that  decays  smoothly 

from  1  to  0  on  the  interval  [0, 1]  so  that  P^^^(O)  —  P^^^(l)  =  0  for  m  =  1, 2, 3,  and  4  (see  Section  2.5).  The 
function  /  of  (4.1)  obviously  has  four  continuous  derivatives  everywhere  with  respect  to  all  its  arguments. 
Notice,  the  temporal  behavior  of  this  source  function  has  been  purposely  chosen  sufficiently  complex;  the 
frequency  of  the  magnitude  oscillations  and  that  associated  with  the  motion  are  incommensurable. 

Equation  (2.15)  driven  by  the  RHS  (4.1)  was  integrated  on  the  fine  grid  of  dimension  512  x  1024  till 
t  =  50  ■  d/c  using  the  lacunae-based  algorithm  of  Section  2  implemented  with  the  fourth-order  scheme  (2.23). 
We  have  chosen  here  a  shorter  time  interval  compared  to  those  we  have  used  for  previous  demonstrations 
(Sections  2.5  and  4.1)  so  that  not  to  make  the  computation  of  the  reference  solution  excessively  expensive. 
This  interval  is  still  quite  sufficient  for  experimentally  judging  the  convergence,  see  Figure  4.3  below.  Having 

all  numerical  examples  we  consider  only  straightforward  motion  because  its  direction  has  to  be  aligned  with  the  z 
direction  of  the  cylindrical  coordinate  system;  otherwise,  the  symmetry  will  be  lost.  In  general,  this  is  not  a  limitation. 
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computed  the  fine- grid  reference  solution,  we  then  integrated  the  same  problem  on  a  collection  of  coarser 
grids  using  the  ABCs  and  compared  the  results  with  this  reference  solution.  In  so  doing,  we  have  employed 
only  the  tw'o  node-centered  schemes:  The  second-order  scheme  (2.22)  and  the  fourth-order  scheme  (2.23). 
The  reason  is  that  when  both  the  fine-grid  solution  and  the  coarse-grid  solution  are  calculated  using  a  node- 
centered  scheme,  it  is  very  easy  to  compare  them  point-wise  (e.g.,  taking  every  other,  every  fourth,  etc., 
node  of  the  fine  grid).  In  contradistinction  to  that,  if  we  were  to  calculate  a  coarser-grid  solution  using 
the  cell-centered  scheme  (2.21),  then  to  compare  it  against  the  reference  solution  we  would  have  had  to  use 
interpolation  on  the  grid.  This  has  a  potential  of  contaminating  the  results  because  of  the  interpolation 
error,  therefore  we  did  not  perform  the  aforementioned  comparison  for  scheme  (2.21). 


Lacunae-based  ABCs  for  (2,2)  node-centered  scheme  Lacunae-based  ABCs  for  (2,4)  node-centered  scheme 


(a)  The  second-order  scheme  (2.22)  (b)  The  fourth-order  scheme  (2.23) 

Fig.  4.3.  Convergence  of  the  solution  of  equation  (2.15),  (4-^)  obtained  with  the  ABCs  to  the  fine-grid  reference  solution. 


Lacunae-based  (2,2)  node-centered  scheme  Lacunae-based  (2,4)  node-centered  scheme 


(a)  The  second-order  scheme  (2.22) 


(b)  The  fourth-order  scheme  (2.23) 


Fig.  4.4.  Same  as  Figure  4.3,  but  the  solution  of  (2.15),  (4.1)  is  obtained  with  the  original  lacunae-based  algorithm. 
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When  computing  the  solution  of  equation  (2.15)  driven  by  non-uniformly  moving  source  (4.1),  we,  of 
course,  employ  the  computational  domain  S{t)  that  traces  the  motion  of  the  source.  As  such,  the  ABCs  are 
set  on  an  artificial  boundary  that  performs  accelerated  motion.  On  Figure  4.3,  we  compare  the  solutions 
obtained  with  the  help  of  the  lacunae-based  ABCs  with  the  fine-grid  reference  solution.  Figure  4.3  shows 
the  dependency  of  the  numerical  error  on  dimensionless  time  for  different  computational  variants.  As  before, 
we  use  the  sequence  of  three  grids:  64  x  128,  128  x  256,  and  256  x  512,  to  demonstrate  the  convergence. 
Figure  4.3(a)  clearly  indicates  the  second-order  convergence  of  scheme  (2.22).  For  scheme  (2.23),  we  observe 
the  fourth-order  convergence  on  Figure  4.3(b). 

To  assess  the  performance  of  the  boundary  conditions,  we  have  computed  the  same  solution  on  the  same 
collection  of  coarser  grids  (64  x  128,  128  x  256,  and  256  x  512)  with  the  same  schemes  (2.22)  and  (2.23),  but 
with  no  ABCs,  rather  using  the  original  lacunae-based  algorithm  of  Section  2  (as  we  did  when  computing 
the  reference  solution).  On  Figure  4.4,  we  compare  the  results  with  the  exact  solution.  As  expected,  scheme 
(2.22)  converges  uniformly  in  time  with  the  second  order  (see  Figure  4.4(a)),  and  scheme  (2.23)  —  with  the 
fourth  order  (see  Figure  4.4(b)). 

Comparing  Figure  4.3(a)  against  Figure  4.4(a)  we  conclude  that  for  the  second-order  scheme  (2.22), 
the  introduction  of  the  ABCs  again  gives  rise  to  no  reflection  back  into  the  computational  domain  (i.e.,  no 
reflection  beyond  the  level  of  the  truncation  error  in  the  interior,  cf.  Section  4,1).  As  concerns  the  fourth- 
order  scheme  (2.23),  one  can  still  notice  slight  differences  between  the  respective  curves  on  Figures  4.3(b)  and 
4.4(b).  The  difference  is  most  visible  for  the  finest  grid  256  x  512,  less  visible  for  the  medium  grid  128  x  256, 
and  non-existent  for  the  coarsest  grid  64  x  128.  This  indicates  that  a  small  amount  of  reflections  due  to  the 
ABCs  may  be  present  in  the  solution,  although  the  actual  elevation  of  the  error  on  Figure  4.4(b)  compared  to 
4.4(a)  is  so  low  that  we  can  regard  these  reflections  negligible  anyway.  Nonetheless,  the  discrepancy  between 
the  corresponding  curves  needs  to  be  accounted  for.  We  attribute  it  to  the  higher  sensitivity  of  the  fourth- 
order  algorithm  to  the  quality  of  the  discrete  lacunae.  This  phenomenon  is  commented  on  in  Section  4.4;  it  is 
not  of  the  fundamental  nature,  the  quality  of  the  lacunae  can  rather  be  controlled  by  appropriately  choosing 
the  parameters  of  the  numerical  procedure,  more  precisely,  the  shape  and  smoothness  of  the  multiplier  fi. 

Summarizing  for  the  case  of  accelerated  motion,  we  see  that  the  solution  of  non-deteriorating  quality 
on  long  time  intervals  can  still  be  successfully  computed  using  the  lacunae-based  ABCs.  To  the  best  of 
our  knowledge,  no  other  ABCs’  methodology  available  in  the  literature  can  handle  artificial  boundaries  of 
domains  that  move  with  acceleration  while  always  keeping  the  reflections  on  or  below  the  level  of  truncation 
error  that  pertains  to  a  given  interior  discretization. 

4.3.  Variable  Speed  of  Sound.  For  the  last  set  of  numerical  experiments  that  we  present  in  the 
paper,  we  wanted  to  pick  up  a  case  that  would  supposedly  be  more  prone  to  the  buildup  of  numerical  error 
inside  the  computational  domain  S{t).  At  the  same  time,  we  wanted  to  keep  the  computational  setup  “in 
line”  with  the  previous  experiments,  see  Sections  2.5,  4.1,  and  4.2.  In  this  connection,  we  notice  that  all  the 
schemes  that  we  have  been  using  for  numerical  demonstrations  previously  were  of  explicit  central-difference 
type.  Consequently,  the  associated  discretization  error  was  of  primarily  dispersive  nature.  In  the  example  of 
the  current  section,  we  wdll  artificially  increase  the  numerical  dispersion  inside  the  computational  domain  S{t) 
and  experimentally  assess  the  performance  of  the  combined  methodology  (interior  scheme  and  the  ABCs)  for 
this  case.  We  should  note,  however,  that  the  use  of  the  ABCs  is  certainly  not  limited  to  the  aforementioned 
class  of  the  schemes. 

It  is  known  that  numerical  dispersion  for  central-difference  schemes  is  more  visible  for  more  “suboptimal” 


29 


Courant  numbers.  In  other  words,  the  further  below  the  stability  limit  the  Courant  number  is,  the  more 
dispersive  the  numerical  waves  become.  In  particular,  it  is  easy  to  see  that  the  one-dimensional  second- 
order  scheme  (3.2b)  is  exact,  and  simply  reduces  to  pure  propagation  along  the  characteristics,  when  the 
Courant  number  ^  is  equal  to  1.  Reducing  this  number  will  introduce  dispersion  of  numerical  waves.  (Of 
course,  the  convergence  of  the  scheme  still  implies  that  the  phase  shift  for  every  given  frequency  will  become 
smaller  as  the  grid  sizes  becomes  smaller.)  The  analysis  of  the  one-dimensional  case  also  indicates  that  in 
multi-dimensional  settings  numerical  dispersion  is  unavoidable.  This  is  easy  to  understand  already  from 
the  following  qualitative  consideration:  To  guarantee  stability  for  all  the  waves  propagating  at  an  angle 
with  respect  to  the  grid  lines  one  has  to  choose  a  smaller  Courant  number,  which  will  necessarily  appear 
suboptimal  for  those  weaves  that  propagate  along  the  grid  lines. 

As  has  been  mentioned,  our  intention  now  is  to  increase  the  numerical  dispersion  inside  the  computational 
domain  and  subsequently  test  the  performance  of  the  combined  algorithm.  To  do  that,  we  continuously 
reduce  the  speed  of  sound  c  from  the  peripheral  areas  of  S{t)  to  its  center.  As  stability  across  the  entire 
domain  will  still  be  limited  by  the  maximum  speed  of  sound,  the  corresponding  Courant  number  near  the 
center  will  be  suboptimal.  This  will  imply  higher  levels  of  dispersion  closer  to  the  domain  center.  This  will 
also  mean  that  any  perturbation  that  originates  in  S{t)  will  stay  inside  the  domain  longer  compared  to  the 
previously  analyzed  cases  of  constant  c.  The  explanation  is  obvious  —  the  interior  speed  of  propagation  is 
low^er.  Consequently,  we  may  expect  that  every  particular  wave  will  accumulate  more  error  before  it  leaves 
the  domain  S{t). 

We  emphasize  that  accurate  quantification  of  the  aforementioned  phenomena  is  not  of  the  central  interest 
for  discussion  in  the  current  paper.  Therefore,  we  do  not  attempt  to  quantify  the  foregoing  considerations, 
especially  as  it  may  appear  technically  difiicult  in  any  non-trivial  setting.  However,  even  on  the  level  of  qual¬ 
itative  understanding  of  the  mechanisms  of  numerical  dispersion,  it  is  certainly  of  interest  to  experimentally 
assess  the  performance  of  the  scheme  with  the  lacunae-based  ABCs  for  the  case  of  variable  c. 

For  our  actual  computations,  we  have  chosen  the  following  law  of  variation  for  the  speed  of  sound  inside 
the  computational  domain  S{t): 


(4.2) 


w^here  c  is  the  original  constant  speed  of  sound  in  the  far  field,  d  denotes  the  diameter  of  the  computational 
domain  as  before,  and  the  polynomial  P(‘)?  well  as  the  quantities  r  and  k,  have  been  introduced  in 
Section  2.5.  The  constant  Pq  in  expression  (4.2)  determines  the  extent  of  reduction  in  the  speed  of  sound  at 
the  center  of  5(t);  and  we  have  tried  two  specific  values:  Pq  =  0.9  and  Pq  =  0.99,  in  our  simulations.  The 
solution  that  we  have  been  computing  in  the  case  of  variable  speed  of  sound  is  the  same  traveling-source 
solution  (2.19b),  (2.20)  that  we  analyzed  before.  However,  instead  of  equation  (2.15)  we  are  now  solving 


df- 


\r  dr  \  dr  J  dz^-  ) 


f{r,z,t)  ,  t>0, 


(4.3) 


w^here  is  defined  by  (4.2).  Substituting  </?(r,  z,  t)  of  (2.20)  into  the  left  hand  side  of  equation  (4.3)  we  obtain 
the  source  term  f{r,z,t)  that  will  obviously  differ  from  that  used  previously  in  Sections  2.5  and  4.1.  This 
new’  source  term  is  still  compactly  supported  on  the  domain  S{t)  C  for  ail  times,  and  it  now  drives  the 
solution  (2.20)  on  the  entire  space.  As  concerns  the  methodology  for  setting  the  ABCs,  it  remains  exactly 
the  same  as  before.  Indeed,  we  point  out  that  the  variation  of  the  speed  of  sound  pertains  to  the  original 
interior  problem  only.  And  the  auxiliary  problem  that  we  solve  for  the  purpose  of  setting  the  ABCs  is  by 
definition  formulated  with  the  constant  speed  of  sound  throughout  its  entire  domain. 
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Lacunae-based  ABCs  for  (2,2)  cell-centered  scheme  Lacunae-based  ABCs  for  (2,2)  node-centered  scheme 


Variable  speed  of  sound  inside  the  computational  domain  Variable  speed  of  sound  inside  the  computational  domain 


dimensionless  time  dimensionless  time 

(a)  The  second-order  scheme  (2.21)  (b)  The  second-order  scheme  (2.22) 

Fig.  4.5.  Convergence  of  the  solution  of  equation  (4.3),  (4.2)  with  the  ABCs  to  the  exact  solution  (2.20)  for  Po  =  0.9. 


Lacunae-based  ABCs  for  (2,4)  node-centered  scheme 


Variable  speed  of  sound  inside  the  computational  domain 


dimensionless  time 


Fig.  4.6.  Same  as  Fig.  4-5  for  the  fourth-order  scheme  (2.. 


In  Figures  4.5  and  4.6  we  present  the  results 
of  the  convergence  study  for  the  case  Pq  =  0.9 
(see  formula  (4.2))  on  the  same  sequence  of  three 
grids  that  we  have  used  previously:  64  x  128, 
128  X  256,  and  256  x  512.  From  Figures  4.5(a) 
and  (b)  we  conclude,  as  before,  that  the  algorithm 
converges  with  the  second  order  for  both  schemes 
(2.21)  and  (2.22);  and  on  Figure  4.6  we  again 
observe  the  fourth-order  convergence  of  scheme 
(2.23).  The  convergence  obviously  does  not  dete¬ 
riorate  as  the  time  elapses,  at  least  till  dimension¬ 
less  time  reaches  the  moment  200  •  djc  when  we 
stop  the  computation.  This  shows  that  similarly 
to  the  previous  cases  of  constant  c,  the  proposed 
numerical  procedure  in  the  case  of  variable  speed 
of  sound  is  still  capable  of  providing  the  solution  of 


non-decreasing  quality.  Note,  however,  that  as  the  original  lacunae-based  algorithm  cannot  be  applied  to  the 
equation  with  variable  c,  we  cannot  directly  compare  in  this  case  numerical  results  obtained  with  the  ABCs 
against  those  obtained  with  no  ABCs  as  we  did  before  (see,  e.g.,  comparison  of  the  results  on  Figure  4.3 
with  those  on  Figure  4.4  in  Section  4.2).  As  such,  in  making  a  conclusion  that  the  ABCs  in  this  case  perform 
practically  as  well  as  they  did  in  the  previous  cases,  we  rely  on  the  foregoing  experimental  observation  of 
temporally  uniform  convergence,  as  well  as  on  the  fact  that  the  actual  error  levels  on  Figures  4.5  and  4.6  are 
only  slightly  higher  than  the  respective  levels  on  Figures  4.1  and  4.2.  This  is  expected,  because  the  results 
on  Figures  4.1  and  4.2  correspond  to  numerically  reproducing  the  exact  same  solution  (p{r,z,t)  of  (2.20) 
using  the  ABCs  but  applying  them  to  the  original  constant-coefficient  wave  equation  in  the  interior. 

Similar  conclusions  as  to  the  convergence  and  quality  of  the  numerical  solution  can  be  drawn  for  the  case 


pQ  =  0.99  from  looking  at  Figures  4.7(a)  and  (b)  and  Figure  4.8.  Prom  the  qualitative  considerations  above 
it  follows  that  this  case  is  supposed  to  be  “tougher”  to  compute,  because  the  reduction  in  the  speed  of  sound 
is  more  significant.  In  practice,  this  is  manifested  by  noticeably  more  oscillatory  error  profiles,  although  we 
still  clearly  see  that  there  is  no  deterioration  of  the  solution  in  the  long  run.  Besides,  the  actual  levels  of  the 
error  are  somewhat  higher  compared  to  the  corresponding  curves  on  Figures  4.5(a)  and  (b)  and  Figure  4.6. 
This  is  also  expected  because  the  numerical  dispersion  inside  S{t)  is  supposed  to  be  higher  as  well. 


Lacunae-based  ABCs  for  (2,2)  cell-centered  scheme  Lacunae-based  ABCs  for  (2,2)  node-centered  scheme 
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(a)  The  second-order  scheme  (2.21)  (b)  The  second-order  scheme  (2.22) 

Fig.  4.7.  Convergence  of  the  solution  of  equation  (4-3),  (4'^)  to  the  exact  solution  (2.20)  for  Pq  =0.99. 


Lacunae-based  ABCs  for  (2,4)  node-centered  scheme 
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4.4.  Implementation  Notes.  The  forego¬ 
ing  algorithm  of  lacunae-based  ABCs  has  several 
parameters  that  need  to  be  tuned  appropriately  so 
that  to  obtain  the  best  possible  results.  Most  of 
the  flexibility  associated  with  the  algorithm  resides 
in  constructing  the  multipliers  and  artificial  near¬ 
boundary  sources  needed  for  computing  the  ABCs 
(Section  3),  as  well  as  in  choosing  the  parameters 
of  the  lacunae-based  integration  (Section  2).  We 
have  not  yet  conducted  a  comprehensive  study  of 
how  the  corresponding  parameters  affect  the  nu¬ 
merical  procedure  and  as  such,  will  only  outline 
here  some  general  trends. 


As  has  been  mentioned  before,  the  multiplier 
Fig.  4.8.  Same  as  Fig.  4.7  for  the  fourth-order  scheme  (2.23).  smooth  in  the  transition  region,  see  Fig¬ 

ures  3,2  and  3.3.  Otherwise,  lacunae  of  the  contin¬ 
uous  solution  will  not  be  reproduced  sufficiently  accurately  in  the  discrete  solution  of  the  auxiliary  problem 
(essentially  because  the  scheme  will  lose  consistency,  see  discussion  in  the  end  of  Section  3.2).  In  most 
of  our  computations  we  have  used  an  algebraic  polynomial  function  with  four  continuous  derivatives  for 
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multiplier,  and  the  extent  of  the  transition  region  was  about  ten  grid  cells.  This  has  always  been  sufficient 
for  the  second-order  schemes.  In  other  words,  we  could  always  obtain  the  temporally  uniform  second-order 
convergence  with  these  settings,  although  we  did  not  check,  for  example,  whether  or  not  it  was  possible 
to  further  reduce  the  extent  of  the  transition  region.  As  concerns  the  fourth-order  convergence,  we  did  see 
situations  when  the  previous  settings  turned  out  somewhat  insufficient  (e.g.,  in  Section  4.2).  This  occurred 
mostly  when  going  from  the  medium  grid  128  x  256  to  the  finest  grid  256  x  512.  To  maintain  the  convergence 
rate  in  this  case,  we  had  to  use  a  wider  transition  region  (fifteen  cells)  and/or  a  smoother  multiplier  (five 
continuous  derivatives).  This  indicates  that  in  general  the  algorithm  of  lacunae-based  ABCs  is  sensitive  to 
the  smoothness  of  the  multiplier,  as  it  is  supposed  to  be.  However,  this  sensitivity  does  not  actually  manifest 
itself  before  the  error  reaches  sufficiently  low  levels.  As  such,  in  practical  computing  one  will  most  likely  be 
able  to  use  rather  narrow  transition  regions,  as  well  as  multipliers  with  limited  smoothness. 

As  concerns  choosing  the  parameters  of  the  lacunae-based  integration  (see  Section  2),  there  is  at  least  one 
important  observation  that  has  been  made  experimentally.  In  theory,  the  contribution  of  a  given  fragment 
of  the  RHS  can  be  subtracted  from  the  overall  solution  as  soon  as  the  time  interval  Tint  has  elapsed  since 
its  inception.  In  practice,  it  has  been  found  useful  to  introduce  the  so-called  aft  front  time  gap  <5,  i.e.,  allow 
for  a  little  extra  time  for  the  waves  to  propagate  outward.  This  implies  choosing  a  somewhat  larger  value 
for  the  period  Z  compared  to  the  necessary  minimum  given  by  (2.6)  so  that  by  the  time  of  subtraction, 
which  is  Tint  +  S  >  Tint,  the  reflected  waves  will  not  have  started  re-entering  the  domain  S{t)  yet.  Moreover, 
we  can  choose  to  introduce  the  actual  front  time  gap  as  well,  i.e.,  increase  Z  even  further,  so  that  by  the 
time  of  subtraction  of  a  given  contribution  the  corresponding  reflected  waves  still  be  at  a  (small)  distance 
from  S{t)  rather  than  right  next  to  its  boundary.  Experimentally,  we  have  found  that  the  aft  front  time  gap 
affects  the  quality  of  the  solution  stronger  than  the  actual  front  time  gap.  However,  the  quantity  S  in  all 
our  simulations  was  sufficiently  small  anyway,  about  4%  of  Tint,  and  most  likely  it  could  be  reduced  even 
further. 


5.  Conclusions.  We  have  constructed  and  tested  the  algorithm  for  setting  highly-accurate  global  arti¬ 
ficial  boundary  conditions  in  the  problems  of  time-dependent  wave  propagation.  The  key  building  block  of 
the  new  ABCs  is  a  special  non-deteriorating  numerical  procedure  that  has  been  developed  previously  for  the 
long-term  integration  of  wave-radiation  problems.  The  latter  procedure  is  based  on  the  presence  of  lacunae 
(aft  fronts  of  the  waves)  in  the  three-dimensional  wave-type  solutions.  The  resulting  lacunae-based  ABCs 
are  obtained  directly  for  the  discrete  formulation  of  the  problem  and  can  complement  any  consistent  and 
stable  finite-difference  scheme.  In  so  doing,  neither  a  rational  approximation  of  non-reflecting  kernels,  nor 
discretization  of  the  continuous  boundary  conditions  is  required.  The  extent  of  temporal  nonlocality  of  the 
new  ABCs  appears  fixed  and  limited,  and  this  is  not  a  result  of  any  approximation  but  rather  a  direct  con¬ 
sequence  of  the  fundamental  properties  of  the  solution.  The  proposed  ABCs  can  handle  artificial  boundaries 
of  irregular  shape  on  regular  grids  with  no  fitting/adaptation  needed.  Besides,  they  possess  a  unique  capa¬ 
bility  of  being  able  to  handle  boundaries  of  moving  computational  domains,  including  the  case  of  accelerated 
motion.  We  have  conducted  a  series  of  numerical  experiments  that  would  corroborate  the  theoretical  design 
properties  of  the  algorithm.  The  experiments  included  computation  of  unsteady  wave-radiation  solutions 
over  long  time  intervals.  In  all  our  experiments  the  ABCs  could  always  keep  the  level  of  reflections  from  the 
artificial  boundary  on  or  below  the  level  of  truncation  error  for  the  interior  discretization  for  as  long  as  the 
computation  was  run.  Besides  the  classical  wave  equation  that  we  have  analyzed  in  the  current  paper,  the 
proposed  technique  may  find  applications  in  computational  acoustics  and  computational  electromagnetics. 
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